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1  Introduction 


Friction  and  rubbing  of  materials  are  among  the  most  common  phenomena  in  mechanics, 
present  whenever  two  solid  bodies  come  into  contact.  They  are  responsible  for  a  variety  of 
occurrences  in  everyday  life.  Some  of  them,  such  as  tire  traction,  are  very  useful;  others,  like 
violin  music,  are  asthetic  and  pleasing,  and  many  others,  such  as  noises,  vibrations,  and  wear, 
are  extremely  unpleasant  and  dilatorius  to  mechanical  systems.  This  common  occurrence 
of  friction  and  the  diversity  of  its  effects  underscore  the  extreme  importance  of  a  deep 
understanding,  modeling,  and  the  control  of  friction  phenomena.  It  is  well  known,  however, 
that  the  phenomena  of  contact  and  friction  of  solid  bodies  are  among  the  most  complex 
and  difficult  to  model  of  all  mechanical  events,  primarily  due  to  the  complex  structure  of 
engineering  surfaces,  the  severe  elasto-plastic  deformation,  damage,  heat  generation,  atomic- 
range  interactions  that  take  place  on  typical  contact  surfaces,  the  presence  of  contaminants, 
lubrication,  and  even  chemical  reactions  on  these  contact  surfaces. 

Efforts  toward  an  understanding  of  frictional  phenomena  and  of  modeling  friction  began 
with  the  historical  works  of  Amontons  [2]  and  Coulomb  [32]  over  two  centuries  ago.  Since 
then,  an  extensive  body  of  experimental  and  theoretical  work  has  accumulated  on  general 
tribology  and  a  good  empirical  understanding  of  the  subject  exists  today.  However,  the 
progress  in  formulating  a  theoretical  background  and  models  of  frictional  interfaces  have 
been  much  slower  to  evolve  than  experimental  investigations.  Although  considerable  progress 
in  this  direction  has  been  made  in  recent  years,  there  are  still  several  issues  that  need 
to  be  resolved  in  order  to  model  friction  and  predict  frictional  phenomena  with  practical 
reliability.  One  of  the  most  difficult  problems  encountered  is  the  estimation  of  material 
constants  occurring  in  the  new  constitutive  models  of  frictional  interfaces.  These  difficulties 
reflect  an  urgent  need  for  constructing  new  constitutive  models  of  contact  and  friction  and 
for  estimating  the  necessary  material  coefficients. 

Presently  there  are  two  basic  approaches  for  the  development  of  mechanical  constitutive 
models  of  friction  and  two  resulting  types  of  frictional  interface  models  which  appear  to  be 
the  most  promising.  There  are; 

1.  phenomenological  models  based  primarily  on  experimental  observations,  and 

2.  a.sperity-bcised  models,  formulated  via  a  theoretical  analysis  and  statistical  homoge¬ 
nization  of  the  microscale  deformation  of  surface  cisperities  in  contact  with  an  opposing 
surface. 

Unfortunately,  to  date  none  of  these  approaches  has  produced  completely  satisfactory 
results.  It  is  well  known  that  experimental  results  depend  strongly  on  the  characteristics  of 
the  test  apparati,  so  the  results  of  different  tests  on  the  same  sample  can  be  considerably 
scattered.  Moreover,  these  macroscopic  experimental  measurements  do  not  provide  sufficient 
insight  into  the  nature  of  the  phenomena  occurring  on  the  contacting  surfaces,  severity  of 
the  deformation,  propagation  of  damage,  etc.  Such  an  insight  can  be  provided  by  asperity- 
based  models,  which  are  based  on  the  microscale  analysis  of  deformation  and  the  relative 
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sliding  of  surface  asperities.  However,  predictions  made  using  classical  or  existing  asperity- 
based  models  are  not  generally  applicable  to  the  environments  normally  met  in  engineering 
applications.  The  main  reason  is  that  these  classical  models  are  based  on  analytical,  closed- 
form  solutions  of  the  deformation  of  a  surface  asperity,  which  require  gross  simplifications  of 
geometry  of  the  asperity  and  of  the  constitutive  models  of  the  contacting  materials  (elastic, 
plastic,  or  at  most  elasto-plastic). 

1.1  Summary  of  Objectives  of  the  Project 

In  this  project,  a  new  approach  for  constructing  constitutive  models  of  frictional  interfaces  is 
under  development,  which  promises  to  provide  a  realistic  link  between  microscale  phenomena 
occurring  on  contacting  surfaces  and  macroscale  phenomenological  models  of  the  interface. 
This  approach  involves  the  use  of  special  finite  element  methods  in  the  modeling  of  complex 
deformations  of  asperities  of  arbitrary  shape,  with  realistic  nonlinear  constitutive  models  of 
the  contacting  materials.  The  technical  approach  for  the  evaluation  of  the  microasperity- 
based  models  of  contact  and  friction  consists  of  two  stages: 

1.  Apply  the  finite  element  technique  to  the  analyze  the  nonlinear  mechanical  responses 
of  surface  asperities  of  different  heights,  shapes,  and  with  general  viscoelastoplastic 
material  properties  and  with  material  damage. 

2.  Apply  statistical  homogenization  techniques  to  evaluate  macroscopic,  phenomenologi¬ 
cal  constitutive  models  of  the  interface. 

The  approach  being  developed  here  should  provide  a  means  for  generating  a  variety  of 
new  and  useful  models  of  frictional  interfaces.  Depending  on  the  selected  level  of  complexity 
of  the  model  of  the  asperity,  a  viscoeiaistoplastic,  hyperelastic,  or  brittle  material  can  be 
considered,  evolution  of  damage  of  the  surface  can  be  modeled,  and  effects  of  lubrication 
and  surface  contamination  can  be  taken  into  account. 


1.2  Research  Progress 

The  results  of  the  first  year  of  effort  on  this  project  included  a  complete  formulation  of  the 
theory  and  a  corresponding  software  package  for  statistical  homogenization  of  the  interface 
parameters. 

in  the  second  year  of  the  project  the  effort  has  focused  on  the  development  of  a  finite 
element  code  for  modeling  the  nonelastic  surface  asperities,  as  well  as  on  the  design  and 
performance  of  the  verification  experiments,  followed  by  initial  verifications  of  the  numerical 
models  for  the  surface  asperities.  In  particular,  the  following  tasks  were  completed  during 
the  second  year: 

1.  development  of  a  three-dimensional  adaptive  finite  element  code  for  the  analysis  of 
surface  asperities  in  contact  with  opposing  surfaces.  The  starting  point  for  this  effort 
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was  an  existing  in-house  finite  element  kernel,  which  was  extended  and  customized  to 
satisfy  the  objectives  of  this  project.  Our  effort  during  the  Icist  year  has  been  focused 
on  the  implementation  of  elcistic  and  viscoplastic  three-dimensional  solid  models,  the 
development  of  contact  and  sliding  resistance  algorithms  and  on  extensions  of  the 
graphics,  user  interface  and  adaptive  algorithms  specifically  related  to  this  project. 

2.  Design,  development  and  performance  of  Phase  I  of  the  verification  experiments,  ori¬ 
ented  toward  testing  the  numerical  models  of  nonelastic  surface  asperities  in  contact 
with  a  rigid  flat. 

3.  Verification  of  the  finite  element  asperity  models  by  comparison  with  analytical  and 
experimental  results.  The  numerical  predictions  were  compared  with  existing  analytical 
solutions  for  selected  simplified  cases  (Hertz  problem)  and  with  the  experimental  results 
obtained  for  fully  nonlinear,  elastoplastic  contact  problems. 

4.  Introductory  tests  of  the  complete  homogenization  procedure,  which  are  based  on  a 
combination  of  numerical  models  of  the  asperity  with  statistical  homogenization,  to 
predict  the  response  of  random  surfaces  to  contact  and  friction  loads.  Comparisons 
with  existing  simplified  asperity-based  models  (Greenwood-Williamson  [44]). 

1.3  Outline  of  the  Report 

This  report  presents  the  results  of  the  second  year  of  effort  on  this  project,  as  well  as  a  brief 
compilation  of  the  most  important  results  of  year  I.  In  particular.  Section  2  presents  a  study 
of  statistical  methods  of  homogenization  of  interface  parameters.  Of  particular  interest  are 
such  issues  as  extraction  of  surface  statistics  from  profilometric  data,  calculation  of  an  asper¬ 
ity  distribution  for  random  surfaces,  and  the  practical  calculation  of  expected  macroscopic 
parameters  from  a  microasperity  analysis.  In  Section  3,  a  detailed  formulation  of  the  bound¬ 
ary  value  problem  representing  the  deformation  of  a  surface  asperity  is  developed.  This 
formulation  includes  elastic  and  viscoelastoplastic  material  properties,  damage  modeling, 
a  nonpenetration  condition  on  the  contact  plane,  and  boundary  conditions  resulting  from 
adhesion  forces  and  sliding  resistance  of  the  interface.  Section  4  presents  the  background 
of  -he  adaptive  finite  element  technology  developed  for  the  analysis  of  the  deformation  of 
a  microasperity.  A  general  background  of  the  /ip-adaptive  finite  element  methodology  is 
discussed  in  this  section,  together  with  a  detailed  presentation  of  the  numerical  algorithms 
used  for  the  solution  of  elastic  and  viscoplastic  contact  problems. 

The  above  theoretical  part  of  the  report  is  followed  by  examples  and  tests  of  the  mi¬ 
croasperity  analysis.  In  particular.  Section  5  presents  some  basic  tests  of  numerical  models  of 
viscoplastic  material  behavior.  Then,  in  Section  6,  finite  element  models  of  asperity  response 
are  verified  by  comparison  with  the  Hertz  solution  and  with  experimental  mecisurements  per¬ 
formed  for  custom-shaped  asperities.  This  section  is  followed  by  an  introductory  example  of  a 
complete  homogenization,  and  comparison  is  made  with  the  classical  Greenwood-Williamson 
contact  model.  Finally,  in  Section  8,  a  research  forecast  is  presented  for  the  third  year  of  the 
project. 
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1.4  Presentations  and  Publications 


Our  research  related  to  new  models  of  contact  and  friction  was  presented  at  the  li3th  ASME 
Winter  Annual  Meeting,  Anaheim,  California,  November  8-13,  1992. 

The  following  papers  were  published  or  submitted  for  publication  in  1992: 

1.  Tworzydlo,  W.W.,  Becker,  E.B.,  and  Oden,  J.T.,  “Numerical  Modeling  of  Friction- 
Induced  Vibrations  and  Dynamics  Instabilities”,  in  Friction-Induced  Vibration, 
Chatter,  Squeal,  and  Chaos,  Ibrahim,  R.A.  and  Soom,  A.,  editors,  ASME  Publi¬ 
cations,  New  York,  1992,  pp.  13-32. 
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2  Asperity— Based  Models  of  Contact  and  Friction 


One  of  the  major  missions  in  tribology  is  the  development  of  constitutive  models  of  frictional 
interfaces.  Throughout  the  decades  a  variety  of  approaches  and  types  of  models  have  been 
developed.  They  can  be  classified  into  several  groups,  including; 

•  models  based  on  experimental  observations, 

•  microasperity-based  models, 

•  phenomenological  models  developed  from  basic  principles  of  mechanics,  and 

•  models  of  the  type  related  to  plasticity  theory. 

It  should  be  noted  here  that  this  distinction  is  only  of  a  general  nature  and  most  of  the 
models  presented  in  the  literature  combine,  in  some  sense,  features  of  more  than  one  of  these 
groups. 

In  this  project  we  focus  on  the  development  of  new  asperity-based  models  of  contact  and 
friction.  These  models  are  aimed  at  the  development  of  constitutive  equations  of  frictional 
interfaces  via  the  statistical  homogenization  of  the  deformation  of  surface  asperities  subject 
to  contact  with  an  opposing  surface.  The  advantage  of  the  asperity-based  models  is  that 
they  provide  good  quantitative  insight  into  the  phenomena  occurring  at  the  interface  and 
predict  additional  information  hardly  available  from  the  experiment -based  laws,  such  as  the 
surface  plasticity  indices,  microfracture  indices,  etc. 

The  first  contact  model  that  was  constructed  to  predict  the  true  contact  area  can  be 
found  in  a  paper  by  .Abbott  and  Firestone  [1],  in  which  the  contact  surface  was  simulated  in 
a  network  of  spheres  that  are  truncated  upon  indentation  into  a  hard  flat.  By  knowing  the 
hardness  of  the  softer  of  the  two  materials  in  contact,  an  estimate  of  the  true  contact  area 
could  be  made,  assuming  perfectly  plastic  deformations. 

An  important  advance  in  development  of  asperity-based  models  of  contact  is  represented 
by  the  pioneering  paper  of  Greenwood  and  Williamson  [44],  in  which  the  rough  surfaces  were 
viewed  as  a  randomly  distributed  population  of  elastic  asperities  with  randomly  distributed 
asperity  heights.  Each  asperity  was  assumed  to  be  spherical  and  elastic  and  its  deformation 
properties  governed  by  the  Hertz  solution  for  elastic  contact.  Experimental  evidence  was 
provided  to  support  the  assertion,  now  widely  held  by  tribologists,  that  for  '  ”  nn. 'ly  isotropic 
engineering  surfaces,  a  Gaussian  distribution  of  asperities  heights  generally  exists.  In  such 
models,  there  are  no  microfrictional  effects  on  the  asperities,  such  effects  leading  to  second- 
order  changes  in  contact  pressure,  a  result  established  nearly  two  decades  earlier  by  Mindlin 
[63].  In  a  related  paper.  Greenwood  and  Tripp  [43]  showed  that  contact  of  two  rough  surfaces 
with  Gaussian  distributions  of  asperity  heights  on  which  asperity  contacts  were  misaligned 
was  equivalent  to  a  single  elastic  surface  with  a  Gaussian  distribution  of  asperity  heights 
impending  on  a  rigid  flat.  The  use  of  such  statistical  representations  of  surface  topography 
has  since  become  a  popular  approach  in  modeling  both  elastic  and  inelastic  contact. 
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The  Greenwood-Williamson  model  was  based  on  the  assumption  that  only  the  asperity 
height  was  a  random  variable,  and  that  the  radius  R  of  each  peak  was  constant.  Several 
generalizations  of  such  random  topography  models  appeared  in  the  literature  of  the  1970s. 
The  paper  of  VVhitehouse  and  Archard  [99]  extends  the  random-asperity  models  to  include 
random  heights  and  curvatures,  and  Nayak  [67]  provided  a  general  approach  to  random  sur¬ 
face  modeling  using  notions  of  joint  probability  distribution  functions.  In  this  same  vein,  we 
mention  the  work  of  Bush,  Gibson,  and  Thomas  [23],  who  derived  a  joint  probability  distri¬ 
bution  density  function  for  random  asperity  heights  and  curvatures  of  a  random  population 
of  elliptic  paraboloids  in  elastic  contact  with  a  smooth  rigid  flat. 

Such  random-microtopography  models  that  employ  a  deterministic  function  for  asperity 
peak  shapes  are  called  asperity  models.  One  source  of  possible  inconsistency  in  such  models 
has  to  do  with  the  fact  that  a  Gaussian  distribution  of  asperity  heights  and  curvatures 
for  a  given  asperity  shape  may  lead  to  a  non-Gaussian  cumulative  probability  distribution 
of  the  surface  height,  an  unrealistic  result  for  most  ‘‘engineering  surfaces.”  This  problem 
was  addressed  by  Hisakado  [48]  and  Hisakado  and  Tsukizoe  [49],  by  assuming  a  Gaussian 
PDF  (Probability  Density  Function)  for  surface  heights,  with  a  given  deterministic  asperity 
shape,  and  then  deriving  the  PDF  for  peak  heights.  Hisakado  [48]  assumed  a  paraboloidal 
asperity  shape  and  Hisakado  and  Tsukizoe  [49]  a  conical  shape.  Francis  [41]  points  out  that 
the  Hisakado  models  may  lead  to  unrealistic  PDFs  for  asperity  heights,  since  they  may  be 
strongly  dependent  on  the  asperity  shape  and  may  become  negative  for  paraboloidal  and 
conical  shapes. 

Extensions  of  asperity-based  models  to  microcontact  deformation  laws  involving  elasto- 
plastic  deformations  were  first  contributed  by  Hisakado  [48].  Hailing  and  Nuri  [46]  account 
for  plastic  deformation  of  the  interface  by  assuming  that  a  rough  surface  deforms  elastically 
while  contacting  a  nonlinearly  elastic  flat,  representing  strain-hardening,  with  each  micro¬ 
contact  defined  by  a  fully-plastic  spherical  indentation.  Significant  generalizations  of  these 
types  of  asperity  models  can  be  found  in  the  detailed  studies  of  Francis  [41],  who  intro¬ 
duces  the  notion  of  the  sum  surface,  discussed  later  in  the  present  work.  This  enables  one 
to  model  Gaussian  engineering  surfaces  with  asperity  shapes  that  a  paraboloidal  only  at 
their  vertices,  but  which  have  random  heights  and  curvatures,  using  the  joint  PDF  of  Nayak 
[67].  Moreover,  Francis  [41]  also  takes  into  account  elastic  and  fully  plastic  deformations, 
with  strain-hardening,  using  functions  determined  empirically  from  spherical  indentations 
of  various  metals.  We  also  mention  that  an  extension  of  the  Greenwood-Williamson  model 
of  spherical  asperities  with  Hertzian  elastic  contact,  constant  radii,  and  random  heights  to 
cases  in  which  a  transition  to  perfectly  plastic  deformations  occur  was  recently  proposed  by 
Chang,  Etsion,  and  Bogy  [26-28]. 

We  note  that  most  of  the  references  cited  above  dealt  with  attempts  to  model  either 
contact  without  sliding  motion,  or  purely  static  or  quasi-static  friction  effects. 

The  asperity-based  models  of  frictional  interfaces  are  constructed  in  five  basic  steps: 

1.  Perform  a  statistical  analysis  for  the  surface  profile  (profiles). 
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2.  calculate  the  surface  statistics  (distribution  of  surface  height,  gradient,  and  curvature) 
from  one  or  more  set  of  profile  data, 

3.  calculate,  from  surface  statistics,  the  probability  distribution  and  density  of  surface 
asperities  of  different  heights  and  (possibly)  peak  shapes, 

4.  calculate,  by  analytical  or  numerical  methods,  responses  of  representatives  of  a  family 
of  surface  asperities  of  different  shapes  to  prescribed  load  programs, 

5.  calculate,  from  asperity  data  and  the  probabilistic  distribution  of  asperities,  the  ex¬ 
pected  values  of  the  interface  response  (normal  and  tangential  forces,  damage,  etc.)  to 
prescribed  load  programs.  This  response  characterizes  constitutive  properties  of  the 
interface. 

Several  variations  of  this  basic  scheme  may  be  derived  for  random  and  deterministic  surfaces, 
isotropic  or  anisotropic  finish,  etc.  In  this  case  a  general  classification  of  surfaces  presented 
(after  Nayak  [66])  in  Fig.  2.1  is  helpful. 

For  practical  purposes,  it  is  reasonable  to  consider  the  following  three  classes; 

(i)  Gaussian  isotropic  surfaces, 

(ii)  Gaussian  anisotropic  surfaces,  and 

(iii)  other  surfaces,  in  particular  deterministic  surfaces  obtained  by  special  finishing  tech¬ 
niques. 

The  flowchart  illustrating  the  homogenization  procedure  for  these  three  groups  is  pre¬ 
sented  in  Fig.  2.2. 
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SOLID  SURFACE 


Figure  2.1:  A  general  classification  of  surfaces. 
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Figure  2.2:  Flowchart  for  statistical  homogenization  of  interfaces 


The  details  of  these  procedures  will  be  discussed  later.  Here  it  is  important  to  observe 
that  for  Gaussian  isotropic  surfaces  it  suffices  to  gather  profilometric  data  along  only  one 
profile  on  the  surface  and  to  consider  asperities  of  axisymmetric  peak  shapes.  For  Gaussian 
anisotropic  surfaces,  however,  one  needs  at  least  three  nonparallel  profiles  and  asperities  of 
different  principal  curvatures  and  orientations.  Finally,  for  non-random  surfaces,  a  full  two- 
dimensional  map  ^  =  z{x,  y)  of  the  surface  may  be  needed,  and  asperities  may  have  various 
deterministic  shapes,  depending  on  the  surface  finish. 

2.1  Microstructure  of  the  Frictional  Interface 

We  begin  by  considering  the  contact  of  two  deformable  bodies,  /  and  II,  over  a  nominal 
contact  area  Aq,  as  illustrated  in  Fig.  2. .3.  An  element  of  unit  nominal  contact  area  is  isolated 
for  study,  as  indicated  in  the  figure.  The  average  stress  vector  S  over  the  unit  contact  area 
has  components  of  force  P  and  Q  normal  and  tangential  to  the  unit  area,  respectively.  The 
situation  is  ec[uivalent  to  that  of  two  typical  coupons  of  surface  material,  one  taken  from  the 
material  near  the  contact  surface  of  each  body,  pressed  together  with  a  force  P  normal  to 
the  tangent  plane  at  the  center  of  the  coupon  interface  and  simultaneously  subjected  to  a 
shear  force  Q  tangent  to  the  plane.  The  bulk  deformations  of  bodies  /  and  II  are  ignored, 
our  aim  being  only  to  characterize  the  mechanical  properties  of  the  contact  interface.  The 
nominal  unit  surfaces  in  contact  are,  for  the  present,  assumed  to  be  initially  flat  and  parallel 
to  one  another. 

It  is  standard  practice  to  depict  the  approximate  profile  of  rough  engineering  surfaces 
with  a  profilometer  or  stylus,  drawn  across  the  surface,  which  generally  yields  a  jagged  profile 
with  an  exaggerated  vertical  scale  of  the  type  shown  in  Fig.  2.4(a).  We  consider  two  such 
opposing  surfaces  1  and  2  which  are  to  ultimately  come  in  contact.  Refrence  planes  defining 
the  mean  asperity  height  of  each  surface  profile  are  established,  and  we  characterize  the 
shape  of  each  profile  by  introducing  functions  Zi  and  Z2,  given  the  height  of  asperities  above 
the  respective  reference  planes,  i.e.,  the  functions  2,  =  Zi{x,y),i  —  1,2,  with  [x.  y)  a  point  in 
the  parallel  mean-height  reference  planes,  define  the  profiles  of  the  rough  material  surfaces 
1  and  2,  respectively.  The  distance  h  between  planes  is  the  separation  of  the  surfaces,  and 
the  distance  between  actual  opposing  material  points  is  denoted  s.  Thus,  at  a  point  [x.y) 
on  the  reference  plane,  we  have 

s  =  h-z  (2.1) 

where  r  is  the  sum  surface  (see  Francis  [41]), 

2  =  Cl  +  22 

Francis  ha.s  pointed  out  that,  from  the  fact  that  the  sum  2  of  the  surface  heights  appears 
in  the  geometric  relation  (2.1),  the  situation  is  equivalent  to  that  of  a  single  deformable 
surface  of  height  2  =  ci  +  22  approaching  a  rigid  flat,  as  suggested  in  Fig.  2.4(b). 

Clearly,  the  undeformed  surfaces  overlap  whenever 

s{x,y)  <  0 
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As  the  normal  load  pressing  the  surfaces  together  increases,  the  separation  h  decreases 
and  at  each  minimum  of  the  function  s  a  microcontact  nucleates  and  expands  due  to  local 
deformation  of  the  surfaces. 

It  is  of  importance  to  note  that  arguments  presented  by  Francis  [41]  are  of  a  purely 
geometric  nature.  From  a  mechanical  point  of  view,  two  major  objections  can  be  raised 
here: 

1.  The  sum  of  two  asperities  (say,  spherical)  in  contact  with  a  rigid  flat  is  not  mechanically 
equivilant  to  two  spheres  in  contact — see  Fig.  2.5. 

In  particular,  the  distance-force  curves  P  =  P{a)  for  the  two  models  are  different. 
Moreover,  the  “asperity”  peak  on  the  sum  surface  corresponding  to  two  spheres  is  not 
spherical.  It  can  be  shown,  however,  using  the  Hertz  solution,  that  these  differences 
vanish  when  the  ratio  of  asperity  radius  R  to  the  contact  radius  r  goes  to  infinity 
(relatively  smooth  surfaces  at  moderate  loads). 

2.  In  the  case  of  contact  with  friction,  the  sum  surface  approach  will  not  model  the  friction 
component  due  to  the  interlocking  of  asperities.  Similarly  as  above,  the  importance  of 
this  effect  diminishes  with  increasing  surface  smoothness. 
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Figure  2.4:  (a)  Profiles  of  opposing  rough  surfaces,  and  (b)  microtopography  of  surfaces  and 
the  equivalent  sum  surface  2. 


1.3 


In  view  of  these  remarks,  the  sum  surface  approach  seems  to  be  correct  and  justified 
for  typical  engineering  surface  finishes  at  moderate  loads.  Note  that  this  condition  is  also 
required  for  the  satisfaction  of  the  assumption  that  separate  microcontacts  do  not  interact 
mechanically  and  that  contacts  do  not  merge. 

It  is  well  known  in  tribology  that  techniques  used  to  produce  engineering  surfaces  usu¬ 
ally  produce  a  Gaussian  distribution  of  the  surface  heights  2,.  Moreover,  the  sum  r  of  two 
Gaussian  surfaces  is  also  Gaussian;  indeed,  Tallian  [89]  points  out  that  if  zi  and  Z2  are  not 
exactly  Gaussian,  their  sum  surface  will  be  closer  to  Gaussian  than  either  surface.  If  the 
shape  of  an  asperity  is  assumed  to  be  paraboloidal,  as  have  been  done  by  several  authors, 
then  the  peak  heights  and  curvatures  are  correlated  random  variables,  with  the  result  that  a 
Gaussian  distriution  of  heights  and  curvatures  may  lead  to  a  cumulative  probability  distri¬ 
bution  of  surface  heights  which  is  non-Gaussian.  This  issue  has  been  studied  by  Hisakado 
[48],  Hisakado  and  Tsukizoe  [49],  and  by  Francis  [41],  who  assert  that  if  the  peak  shape  is 
paraboloidal  only  at  its  vertex,  then  the  ensemble  of  peaks  can  be  made  to  conform  to  the 
Gaussian  distribution. 

2.2  Statistics  of  a  Random  Surface 

There  are  several  methods  of  homogenization  that  can  be  found  in  the  literature.  Few 
have  been  effectively  used  for  describing  nonlinear  frictional  phenomena.  As  one  possible 
technique  we  describe  a  general  approach  inspired  by  the  works  of  Lonquet-Higgins  [56-58], 
Nayak  [67],  and  Francis  [41];  see  also  Chang,  Etsion.  and  Bogy  [26-28].  To  fix  some  of 
these  ideas,  we  note  that  for  a  given  asperity  profile,  one  defines  the  autocorrelation  function 
C(A'’,  Y)  for  the  random  variable  c  =  c(x,  (/)  (the  surface  height), 

C(A'’.y’)  =  lim  —  /  /  z{x,y)z{x  +  X,y  +  Y)dx  dy 

“^■50  ab  Jo  Jo 
6—00 

and  the  power  spectral  density  P{kj;,ky)  as  its  Fourier  transform, 

P{k.,ky)  =  ^  r  r  C(A,y')e.xp  [-i(Afc,  +  Yk,)]dX  dY 

47r^  */— 00*'— 00 

The  power  spectral  moments  are 

/oo  roo 

/  P(k,,ky)k'Jl  dk.dky  (2.2) 

'OO  »/— 00 

and  the  r.m.s.  roughness  cr  is  the  variance, 

/OO  fOO 

/  P{kx,ky)  dkj.  dky 

•00  J  —00 

A  convenient  representation  of  a  continuous  random  surface  is  of  the  form  [57,67]: 

00 

H  cos{xk^„  +  ykyr^  +  £„)  (2.3) 

n=l 
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where  amplitudes  C„,  wave  numbers  kxn  and  ky„,  and  phase  £„  are  random  variables.  It  is 
assumed  that  there  are  an  infinite  number  of  wave  vectors  in  any  area  dkxdky  and  that  £  has 
a  uniform  probability  density  in  the  range  (0,2t).  The  power  spectral  density  is  related  to 
representation  (2.3)  by 

P(K,k,)dk,dk,  =  Wcl 

the  summation  being  over  all  terms  with  {knx,kny)  lying  in  the  area  dkxdky  around  {kx,ky). 
The  power  spectral  moments  rUp,  can  then  be  expressed  as 

=  (2.4) 

^  n 

Similar  definitions  and  representations  as  above  can  be  introduced  for  arbitrary  surface 
profile  z{s),  s  being  a  parameter  on  the  surface.  Of  particular  interest  are  spectral  moments 
of  a  profile  mo,  m2,  and  m4. 


2.3  Calculation  of  Surface  Statistics  From  Profile  Data 

Information  about  the  statistics  of  a  two-dimensional  surface  can  be  effectively  obtained 
from  profilometric  data  for  one  or  more  profiles  on  the  surface.  This  greatly  simplifies  the 
homogenization  procedure  because  both  experimental  measurements  and  statistical  post¬ 
processing  are  much  easier  for  one-dimensional  profiles.  In  this  section  we  discuss  details  of 
these  computations  for  both  isotropic  and  anisotropic  surfaces. 

2.3.1  Profiles  on  Gaussian  Isotropic  Surfaces 

It  was  shown  by  Lonquet- Higgins  [56,  57]  and  Nayak  [67]  that  for  random  isotropic  surfaces 
the  mean  surface  height  and  non-zero  spectral  moments  are  expressed  in  terms  of  mean 
profile  height  and  profile  spectral  moments: 


^surfare 

—  ^profile 

moo 

=  mo 

ru2o  =  mo2 

=  m2 

3m22  =  mo4  =  m4o 

=  m4 

Therefore,  in  order  to  calculate  surface  statistics  it  suffices  to  perform  mecisurements  for  one 
profile  on  the  surface.  The  spectral  moments  of  the  profile  can  be  calculated  in  several  ways: 

•  from  the  definition  as  moments  of  power  spectral  density 

•  from  statistical  postprocessing  (sampling)  of  profile  data 

•  from  counting  zeros  and  extrema  of  the  profile 
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Calculation  from  definition: 

m.  =  r  ^k)k'dk 
7  —  00 

where  fc  is  a  wave  number  and  $(A:)  is  a  power  spectral  density.  This  particular  method  is 
rather  expensive  because  it  requires  evaluation  of  the  autocorrelation  function  and  the  power 
spectral  density  as  its  Fourier  transform. 

It  is  much  easier  to  calculate  profile  spectral  moments  if  one  reinterprets  them  as  standard 
deviations  cr,  a,  a  of  profile  heights  s,  slopes  i,  and  curvatures  £,  respectively: 

mo  = 
m2  =  (P’ 

1714  =  P 


Assuming  that  the  profile  data  was  sampled  at  n  points  separated  by  the  interval  As  (see 
Fig.  2.6),  the  profile  statistics  can  be  calculated  from  the  following  sampling  formulas  [16]: 

(a)  mean  height,  slope  and  curvature: 


1  " 
rE.-. 


1  ” 

t=i 

1=1 


(b)  variations  of  height,  slope  and  curvature: 


=  Prp-P 

5) 

“  ^  1=1 


The  values  of  first  and  second  derivatives  can  be  calculated  from  a  second  order  approxi¬ 
mation  of  the  profile  shape,  discussed  in  Appendix  A.l  in  reference  [94].  Note  that  the  mean 
slope  and  mean  curvature  of  a  perfect  Gaussian  profile  should  be  zero.  For  real  profiles  they 
may  slightly  differ  from  zero.  Also  note  that  for  the  above  procedure,  shorter  wavelengths 
can  be  automatically  filtered  out  by  appropriate  selection  of  the  sampling  interval  As. 


17 


18 


An  ingenuous  alternative  way  of  calculating  m2  and  was  proposed  by  Lonquet-Higgins 
[57],  see  also  Nayak  [67].  The  densities  of  zeros  and  of  extrema  of  the  profile  are  expressed 
through  spectral  moments  as: 


■^zero 

Dextr 


I 

TT 

TT 


1  a 
TT  cr 

1  a 
TT  a 


By  zero  point  we  mean  here  the  point  where  c(s)  =  z.  The  counting  of  zeros  and  extrema 
can  be  performed  simultaneously  with  the  sampling  procedure  described  above.  Then,  after 
calculation  of  the  variation  cr,  variations  of  slopes  and  curvatures  can  be  obtained  as: 


<T  —  'K  CT  D-xero 
T  =  zaDextr 


In  practice,  both  sampling  and  counting  methods  can  be  easily  implemented  in  the  same 
sampling  program.  Practical  comparisons  of  these  procedures  are  presented  in  our  previous 
report  [94]. 

Another  parameter  necessary  for  the  homogenization  procedure  is  a  density  of  peaks  on 
the  surface,  defined  for  homogenous  surfaces  as: 


lim 

<ix— CO 

dy— 00 


dA 


where  dA  =  dx  dy  is  the  surface  area  and  Np  is  the  number  of  asperity  peaks  within  this 
area.  The  density  of  surface  peaks  can  be  calculated  from  profile  parameters  [67]  as: 


Dr, 


1  rUi 

67r\/5  ^2 


2.3.2  Profiles  on  Gaussian  Anisotropic  Surfaces 

Basic  statistical  information  for  anisotropic  random  surfaces  consists  of  nine  moments  of  the 
power  spectral  density:  moo, ...  ,17140.  However,  since  the  properties  of  the  surface  do  not 
depend  on  the  orientation  of  the  x,  y  axes,  only  certain  invariant  combinations  appear  in  the 
probability  distribution  of  the  surface  statistics  [58,67].  These  invariants  are: 

1.  moo 

2.  mo2  +  m2o 

3.  m2omo2  -  mfi 

4.  m4o  +  2m22  +  ttioa 
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5.  m4oJ^04  —  4mi3m3i  +  3m22 

6.  (jTl40  +  ^^22)  (JT^22  +  ^04)  ■“  (^31  +  ^13)^ 

7.  |m40  (m22mo4  -  mJa)  -  m3i  (m3imo4  -  rniim22)  +  ^22  (^31^13  -  "?|2)  } 

From  three  profiles  in  three  nonparallel  directions  0,,  i  =  1,2,3  nine  parameters  can  be  de¬ 
fined:  mo(,),m2(i),m4(,),f  =  1,2,3.  However,  since  mo(i)  =  mo(2)  =  ”20(3),  then  these  three 
profiles  define  seven  constants — invariants  described  above.  This  means  that  three  nonpar¬ 
allel  profiles  suffice  to  define  surface  statistics  for  Gaussian  anisotropic  surfaces.  (Detailed 
equations  will  not  be  derived  here.) 


2.4  Calculation  of  Asperity  Statistics  From  Surface  Statistics 

The  primary  idea  of  asperity-based  interface  models  is  to  calculate  interface  parameters 
(normal  force,  friction  force,  etc.)  for  a  family  of  asperities  of  certain  deterministic  shapes  and 
to  obtain  expected  values  of  these  parameters  for  the  interface  from  a  statistical  distribution 
of  asperities.  This  requires  the  calculation  of  probability  density  of  surface  asperities.  For 
random  surfaces,  this  probability  density  can  be  expressed  in  terms  of  surface  statistics.  This 
problem  will  be  addressed  in  this  section. 


2.4.1  Asperity  Statistics  for  Gaussian  Isotropic  Surfaces 


For  Gaussian  isotropic  surfaces  two  random  variables  are  assumed  to  govern  the  distribution 
of  asperities:  asperity  peak  height  and  mean  curvature  k.  The  joint  probability  density 
function  of  these  parameters  was  derived  by  Nayak  [66]  and  recast  in  a  different  form  by 
Francis  [41].  Here  we  present  the  formula  due  to  Francis: 


where 


e2< 


a 

nondimensional  peak  height 

T)  =  \/L5  — 

(7 

nondimensional  peak  curvature 

0  = 

ao 

wavelength  spectrum  parameter 

The  wavelength  spectrum  parameter  varies  for  random  surfaces  between  0  and  1;  zero 
corresponds  to  the  widest  wavelength  spectrum  (asperity  heights  and  wave  numbers  are 
not  correlated),  and  one  corresponds  to  the  narrow  spectrum  (longer  asperities  have  bigger 
heights).  Note  that  deterministic  surfaces  may  have  >  1,  see  Section  2.6.1. 
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2.4.2  Asperity  Statistics  for  Gaussian  Anisotropic  Surfaces 

For  Gaussian  anisotropic  surfaces  the  representative  asperities  are  no  longer  of  axisymmetric 
shape.  Instead,  one  should  consider  asperities  with  elliptic  horizontal  cross  sections,  shown 
in  Fig.  2.7b. 

The  peaks  of  these  asperities  can  be  characterized  by  four  parameters: 

1.  2p  —  peak  height 

2.  Ki,K2  —  principal  curvatures 

3.  cr  —  orientation  of  the  main  axis  of  curvature 

Equivalently,  peak  height  Zp  and  three  Cartesian  curvatures  (second  derivatives  of  r)  /c^y, 
Kyy  can  be  used.  The  joint  probability  density  function  based  on  all  these  parameters  should 
be  defined  as  f^piTna  In  principle,  this  function  can  be  defined  from  surface 

statistics,  in  particular  spectral  moments  moo,-..,fn4o  [58,67].  To  the  author’s  knowledge, 
no  such  formula  is  presently  available  in  a  closed  form. 

2.4.3  Asperity  Statistics  for  Deterministic  Surfaces 

Some  special  types  of  finish  may  produce  surfaces  of  non-Gaussian  random  distribution  or 
deterministic  distribution.  For  such  arbitrary  surfaces  the  distribution  of  asperity  peaks 
cannot  be  obtained  from  profile  data  and  need  to  be  calculated  directly  from  a  surface  map. 
In  this  section  we  present  a  simple  sampling  procedure  to  calculate  asperity  statistics  from 
surface  data. 

We  assume  that: 

•  The  surface  is  homogenous. 

•  The  function  z(x,y)  (surface  height)  for  the  surface  is  given.  This  can  be  obtained 
from  two-dimensional  sampling,  holography  or  other  methods. 

•  Asperity  peaks  are  characterized  by  the  peak  height  Zp  and  mean  curvature  k. 

The  probability  density  of  asperity  height  and  curvature  /-^^(^p,  k)  can  be  obtained  from  the 
following  sampling  procedure;  Cover  the  domain  Q  with  a  regular  mesh  of  points  (xi,  yi).  i  = 
l,n  presented  in  Fig.  2.8. 

The  mesh  spacing  h  can  be  defined  so  as  to  filter  out  high-frequency  noise.  The  sampling 
is  performed  by  looping  through  the  points  (x,,yi),i  =  l,n  and  for  each  point: 

1.  Check  if  the  point  is  a  peak  or  near  a  peak  within  resolution  h.  The  peak  is  identified 
as  a  peak  if  its  height  z{Xp,yp)  is  greater  than  all  its  nearest  neighbors  (eight  for  interior 
points).  Alternatively,  more  elaborate  criteria  may  be  used. 
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Figure  2.8:  Sampling  of  arbitrary  surface. 
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2.  If  the  point  is  a  peak,  then  calculate  the  second  derivatives  Zxxi^yy  and  Zxy  Here 
simple  finite  difference  formulas  may  be  used  or  a  generalized  minimization  procedure 
eis  presented  in  Appendix  A. 


The  mean  surface  height  and  standard  deviation  of  the  surface  height  are  calculated  as: 

I  =  -Y- 


The  joint  probability  density  of  asperity  peak  heights  and  curvatures  can  be  calcu¬ 
lated  after  locating  all  the  peaks  by  dividing  the  range  of  peak  heights  and  curvatures 
■^p max]  ^  [^min)  ^max]  into  3,163,  elements  (see  F]g.  2.9). 

Then  for  each  area  element  with  a  center  point  (zpi,Kj)  define 

/zis  (^pii^j)  ~  J) 

Ap 

Here  Np  is  the  total  number  of  peaks  and  np{i,j)  is  the  number  of  peaks  within  the  area 
element  AzpA/c,  identified  by: 


■‘■pmin  "b  l)Az  ^  Zp  <  iAz 

^min  "b  ij  1)A/C  ^  Kp  <  Kmin  "b  JAk 

The  values  above  define  the  discrete  values  of  the  joint  probability  density  function  of  as¬ 
perity  peak  heights  and  curvatures.  This  function  may  then  be  regularized  by  an  application 
of  appropriate  approximation  techniques.  A  similar  procedure  can  be  used  if  one  chooses 
to  characterize  asperity  peaks  with  more  than  two  parmeters,  such  as  Zp,Ki,K2  and  a  for 
anisotropic  surfaces. 


2.5  Calculation  of  Macrocontact  Expectations  of  Interface  Pa¬ 
rameters 

Once  we  know  the  probability  distribution  of  asperity  heights  and  shapes  as  well  as  the 
values  of  any  parameter  X  for  single  asperity,  it  is  possible  to  calculate  the  expected  value 
of  X  for  the  interface.  This  procedure  varies  somewhat  depending  on  the  classification  of 
surface  type.  In  this  section  we  present  a  detailed  procedure  of  the  expectation  calculation 
for  Gaussian  isotropic  surfaces  and  outline  its  extensions  to  other  surface  types. 
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2.5.1  Expectation  Calculation  for  Gaussian  Isotropic  Surfaces 

For  Gaussian  isotropic  surfaces,  the  following  parameters  are  needed  to  calculate  the  expected 
value  of  the  macroscopic  interface  parameter  .V : 

i)  Np  The  number  of  peaks  within  the  contact  area  ^o- 

Probability  density  of  asperity  peaks  of  (nondimensional)  height 
^  and  mean  curvature  i].  By  a  simple  change  of  variables  one 
can  define  f^j,(zp,K). 

iii)  X{zp,  K,a,d)  The  value  of  parameter  .V  for  different  peak  heights  Zp  and  cur¬ 

vatures  K,  subjected  to  a  normal  approach  a  and  sliding  distance 

d. 

iv)  cr,&,iT  Deviations  of  profile  heights,  slopes,  and  curvatures  used  to 

nondimensionalize  peak  heights  and  curvatures. 

The  expected  value  of  X  per  asperity  is  calculated  as 

E(.Y(a,d))  =  j  X{zp,K,a.d)f^^(Zp,K.)dKdzp  (2.5) 

and  the  macrocontact  expectation  of  A*  is 

X{a,d)  =  NpE[X{a.d)) 

Note  that  even  for  Gaussian  isotropic  surfaces  there  exist  asperities  of  nou-axisymmetric 
cross  sections.  However,  due  to  isotropy,  it  suffices  to  consider  only  axisymmetric  represen¬ 
tatives  of  certain  mean  peak  curvature  k.  In  this  project  we  choose  the  asperity  to  be  a 
cosine  hill  defined  in  a  local  coordinate  system  as: 

z(x,  ^)  =  C  cos  kx  cos  ky 

For  this  asperity  the  peak  height  and  mean  curvature  are  defined  as 

=  C 

K  =  Ck'^ 

The  above  shape  is  consistent  with  a  generic  representation  of  a  Gaussian  surface  presented  in 
formula  (2.3).  This  is  different  than  approaches  presented  to  date  in  the  literature,  in  which 
tisperity  peaks  were  usually  assumed  to  be  spherical  or  paraboloidal.  This  was  l)ecause  these 
works  were  based  on  analytical  solutions  for  asperity  deformation,  such  as  Hertz'  solution. 
In  this  work  we  are  modeling  the  asperity  by  the  finite  element  method,  so  it  is  possible  to 
use  the  model  which  does  not  suffer  from  inconsistencies  of  spherical  or  paraboloidal  asperity 
peaks. 

Although  for  Gaussian  isotropic  surfaces  the  probability  density  of  asperity  heights  and 
curvatures  is  analytic,  the  values  of  A^(rp,  k,  . . .)  that  we  obtain  from  finite  element  com¬ 
putations  are  not.  Therefore  the  e.xpected  value  of  .V  must  be  calculated  using  numerical 
quadrature.  This  quadrature  was  implemented  under  the  following  assumptions: 


(a)  The  domain  of  integration  is  truncated  to  the  subregion  [zmin,  •^max]  x  [«min,  ^max], 
where  the  probability  density  is  large  enough  to  effectively  contribute  to  the 
final  integral  E{X).  This  region  is  defined  adaptively  (see  Section  2.6.2). 

(b)  The  parameter  values  X{zp,K, . . .)  are  given  (calculated  by  FEM)  for  certain 
selected  values  of  peak  heights  and  curvatures  in  the  domain  of  integration.  These 
points  are  not  necessarily  regularly  distributed  within  the  domain  of  integration. 

The  numerical  procedure  for  the  calculation  of  the  integral  consists  of  the  following  steps: 

1.  Divide  area  [^mini -s^max]  x  i  ^max]  into  area  elements  Az  x  A/c  (see  Fig.  2.10). 

2.  Calculate  the  integral  by  looping  over  cells  and  applying  numerical  quadrature  (trape¬ 
zoidal,  Simpson,  Gauss,  or  any  other)  according  to  the  formula 

N  r  m 

E{X)  ^  ^  ^  (-^pCor))  ^(0)5  •  •  •)  fin  (-p(o)i  U’o.ArA/v^  (“-6) 

1=1  La=l 

where  i  is  the  number  of  integration  cells,  a  is  the  number  of  quadrature  points  within 
a  cell,  and  Wa  is  the  corresponding  weight  factor. 

Note  that  this  integration  requires  the  value  of  X  at  integration  points  (2:p(a),  «(»))  within 
each  cell.  Since  it  may  be  difficult  to  perform  finite  element  analysis  for  values  of  Zp  and 
K  corresponding  exactly  to  all  the  integration  points,  these  values  are  calculated  from  the 
original  data  points  using  the  error  minimization  procedure  presented  in  our  previous  report 
[94].  The  quadrature  rule  currently  implemented  for  integration  are  the  trapezoidal  and 
four-point  Gauss  rule. 

Note  that  the  above  procedure  introduces  error  due  to  truncation  of  the  integration 
domain  and  due  to  numerical  integration.  This  leads  to  a  rather  unwelcome  result  that  even 
for  constant  A'",  the  calculated  expected  value  £’(A'’)  would  be  different  than  A'.  In  order  to 
compensate  for  this  error,  we  additionally  calculate  the  integral  of  the  probability  density 
(which  should  be  one): 

N  '  m 

I  =  fin  (■^(or)i  ^p(o))  WqAzAk 

i=l  Lo=l 

Then  the  corrected  value  of  expectation  of  A  is  calculated  cis: 

EiX)  =  E(X)/I  (2.7) 

This  procedure  assures  that  for  constant  A  the  expected  value  E{X)  is  equal  to  A'^. 

2.5.2  Expectation  Calculation  for  Random  Anisotropic  Surfaces 

As  mentioned  previously,  for  anisotropic  Gaussian  surfaces  one  has  to  consider  asperities  of 
random  peak  heights  Zp,  principal  curvatures  /ci  and  K2,  and  orientations  of  the  principal 
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axis  a.  The  calculation  of  expected  values  of  the  interface  parameters  is  similar  to  equation 
(2.5); 


Ki ,  /Cj ,  O,  (Z  . 


foo  roo  foo 

I  I  I  I  ^  f^r)irt2Ci  ^2i  dOidKidK2dZp 

J—oo  Jo  Jo  Jo 


(2.8) 


Note  that  presently  there  exist  no  closed  form  solution  for  the  probability  density 
of  asperity  peaks  of  random  heights,  principal  curvatures  and  orientations.  Also  note  that 
in  order  to  span  the  integration  space,  one  would  need  to  obtain  finite  element  solutions  for 
a  large  family  of  asperities,  corresponding  to  various  combinations  of  2p,  /Ci,  /C2,  and  a.  This 
would  be  a  very  expensive  task  computationally,  and  will  not  be  considered  in  this  project. 


2.5.3  Expectation  Calculation  for  Deterministic  Surfaces 

The  calculation  of  expectation  values  E{X)  for  non-random  surfaces  follows  essentially  the 
same  numerical  procedure  as  for  Gaussian  isotropic  or  anisotropic  surfaces.  Depending  on 
the  surface  type,  the  peak  height  Zp,  curvature  k  and  other  parameters  may  be  selected  to 
represent  typical  asperities.  The  joint  probability  density  f^q,„{zp,  k,  . . .)  can  be  obtained 
from  surface  sampling  as  discussed  in  Section  2.4.3. 


2.6  Numerical  Verification  of  Statistical  Postprocessing 

The  homogenization  procedures  discussed  in  the  previous  section  for  Gaussian  isotropic 
surfaces  were  used  as  the  basis  for  the  implementation  of  specialized  software  for  this  purpose. 
Verification  tests  of  this  software  were  presented  in  our  last  yearly  report  [94].  In  this  report, 
we  used  the  above  approach  to  simulate  Greenwood- Williamson  contact  model  (Section  7). 


3  Deformation  Mechanics  of  a  Single  Asperity 

We  now  focus  on  the  analysis  of  a  typical  asperity  in  contact  with  a  rigid  flat.  The  asperity  is 
a  body  of  revolution,  symmetric  about  its  z  =  xa-axis,  and  subjected  to  adhesion  pressures  q 
on  its  exterior  surfaces  that  are  not  in  contact  with  the  rigid  flat,  and  to  contact  pressures  due 
to  its  indentation  into  the  rigid  flat  (see  Fig.  3.12).  The  equations  governing  the  deformation 
of  the  asperity  are  discussed  below. 

3.1  Momentum  and  Geometric  Equations 

The  momentum  equations  for  the  asperity  are; 

=  0  (3.9) 
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where  <Tij  is  the  Cauchy  stress  tensor  at  a  point  x  =  {xi,X2,X3)  €  Jl,  being  the  open 
material  domain  of  the  asperity  and  crijj  is  the  divergence  of  the  stress  cr,j. 

In  rate-dependent  viscoplastic  applications  a  rate  form  of  the  equilibrium  equations  is 
used: 

d.jj  =  0  (3.10) 

where  the  dot  denotes  the  time  derivative. 

Geometric  equations  express  strains  in  terms  of  displacements: 

Strains  can  be  decomposed  into  elastic  and  nonelastic  strains: 

■"‘J  ^IJ  ~  '-tj 

Similarly  as  for  the  momentum  equations,  the  rate  form  of  geometric  equations  will  also 
be  used: 


^ij  "b  ^J.i  ) 


3.2  Constitutive  Equations 

Surface  asperities  for  typical  engineering  surfaces  consist  of  the  same  material  as  the  bulk 
body  with  possible  contaminations  and  structure  change  from  oxidation  and  surface  finish 
processes.  Therefore,  for  general  surfaces  a  variety  of  material  classes  should  be  considered, 
such  as  elastic,  hypoelastic,  elastoplastic,  etc.  In  this  project  two  major  material  classes  are 
considered,  namely: 

•  linearly  elastic  (isotropic  or  anisotropic)  for  some  metal  surfaces,  ceramics,  composites, 
and  hard  rubbers,  and 

•  viscoelastoplastic  models  for  metallic  surfaces  and  modern  ductile  ceramics. 

3.2.1  Linearly  Elastic  Constitutive  Models 

The  general  linearly  elastic  constitutive  relations  are  given  as: 

<^i}  =  EijklSkl 

where  Eijki  are  the  components  of  the  fourth  order  elasticity  tensor.  It  has  up  to  36  inde¬ 
pendent  coefficients  for  general  anisotropic  materials.  However,  for  most  material  classes, 
the  number  of  material  coefficients  is  much  smaller  and,  for  isotropic  materials,  there  are 
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only  two  coefficients,  E  and  i/.  The  specific  forms  of  tensor  E  for  various  materials  are  well 
known  and  will  not  be  presented  here. 


3.2.2  Elasto- Viscoplastic  Constitutive  Model  With  Damage 


We  now  describe  the  Bodner-Partom  constitutive  equations  [10,11]  used  in  the  modeling  of 
viscoelastoplastic  asperities.  The  elastic-viscoplastic  analysis  is  based  on  decomposition  of 
strain  rates 


(3.11) 


where  superscripts  (e)  and  (n)  denote  elastic  and  nonelastic  strain  components,  respec¬ 
tively.  The  constitutive  relations  are 


<T,J  =  -  4?^)  (3.12) 

A  nonelastic  deformation  is  governed  by  the  flow  rule: 

2.  =  gi(<^ij,^k) 

Ui  = 

where  fij,gi  and  h,  are  constitutive  functions,  Zi  are  internal  state  variables,  and  u.’;  are 
damage  variables.  These  functions  and  state  variables  characterize  the  viscoplastic  response 
of  the  material  with  continuum  damage  effects. 

In  the  particular  version  of  the  Bodner-Partom  theory  applied  in  this  work,  the  nonelastic 
flow  rule  is  of  the  form: 


-  As 

tjj  —  AS,j 

where  s,_,  are  the  deviatoric  components  of  a  stress  tensor 

5|j  —  cr,j 

The  current  value  of  parameter  A  is  given  by 


A^  =  exp 


\  3J2  J 


,A  >  0 


where  Ji  is  the  second  invariant  of  a  deviatoric  stress  tensor 
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J2  = 


1 


Dq  is  a  limiting  strain  rate  in  shear,  n  is  a  material  constant  and  2  and  u  are  state  variables 
which  evolve  during  deformation.  In  particular,  z  is  the  hardness  variable,  which  represents 
viscoplastic  hardening  (or  softening)  of  a  material.  The  variable  u  is  the  damage  variable. 
This  variable  represents  weakening  of  the  material  due  to  nucleation  and  propagation  of 
microscopic  voids  and  cracks  in  the  material.  The  micro-cracks  considered  here  are  in  the 
range  of  0.01  mm  in  length.  The  rupture  criterion  is  cj  =  1,  which  corresponds  to  the 
saturation  of  the  material  with  voids.  Alternatively,  a  single  crack  may  grow  to  a  size  on  the 
order  of  1  mm.  In  the  latter  case,  crack  is  too  big  to  be  treated  in  a  continuum  sense,  and 
its  propagation  should  be  followed  using  the  methods  of  fracture  mechanics. 

In  the  framework  of  materials  science  the  value  of  ui  is  usually  interpreted  as  ratio  of  the 
area  of  voids  to  the  total  area  of  a  certain  cross  section  of  a  sample: 


u  — 


^void 


A 


The  state  variables  «  and  u  evolve  according  to  the  specific  equations  of  the  viscoplastic 
theory: 

1.  Evolution  equations  of  hardness  variable 

The  internal  state  variable  2  consists  of  isotropic  and  directional  components, 


2^  -b  2^ 


The  evolution  equation  proposed  for  the  isotropic  hardening  component  [10,11,24,25]  is 


i'(0  =  m.ln  -  /(OIVMO  - 


(t)  —  22 


(.3.13) 


with  the  initial  condition,  2:^(0)  =  zq.  In  the  first  term,  zi  is  the  limiting  (saturation)  value 
of  z^rrii  is  the  hardening  rate,  and  the  plaistic  work  rate  is 

W  =  cr- 

which  is  taken  as  the  measure  of  hardening.  22  is  the  minimum  value  of  2^  at  a  given 
temperature,  and  Ai  and  ri  are  temperature  dependent  materia)  constants.  The  evolution 
form  of  the  directional  hardening  component  (Refs.  [10,11,24,25])  is  defined  as 

where  u,;  are  the  direction  cosines  of  the  current  stress  state. 


Uij(t)  =  <7ij{t)/[(7kiaki]^ 


(.3.14) 
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The  evolution  equation  for  0ij{t)  has  the  same  general  form  as  that  for  isotropic  hardening 
but  has  tensorial  character, 

Aj  =  rn2{z3Uijit)  -  l3ij{t)]Wp{t) 

where 

and 

A,(0)  =  u 

As  in  Eq.  (3.13),  is  the  hardening  rate.  and  r2  are  temperature  dependent  material 
constants. 


2.  Evolution  of  damage 

The  damage  parameter  consists,  in  general,  of  isotropic  and  directional  components, 

U)  =  A 


The  evolution  of  isotropic  damage  proposed  in  reference  [11]  is  of  the  form 


H 


In 


1 


E±i 

P 


(3.15) 


In  the  above  P  and  H  are  material  constants,  Q  is  the  stress  intensity  function,  given  by 


Q  = 


+  B^,  +  Cl* 


where  is  the  maximum  principal  tensile  stress,  1^  is  the  first  stress  invariant  (nonneg¬ 
ative)  and  J2  is  the  previously  introduced  second  invariant  of  deviatoric  stress. 

A,  5,  C,  and  u  are  material  constants.  A,  5,  C  must  satisfy  the  condition 

A  +  B  +  C=^l 


Clearly,  the  actual  proportion  of  these  constants  selects  the  factor  for  stress  state  'vhich  is 
most  important  in  the  development  of  internal  damage. 

The  initial  condition  for  isotropic  damage  is  tt'^(O)  =  0.  In  practical  analyses  the  coef¬ 
ficient  u  is  of  the  order  10  (compare  ref.  [11]).  Thus,  when  SI  (metric)  units  are  used  in 
the  anal3'sis,  the  factor  Q  as  well  as  the  constant  H  reach  extremely  high  values,  beyond  the 
limit  of  real  number  capacity  on  some  computers.  Thus,  for  numerical  analysis,  equation 
(3.15)  was  recast  in  the  equivalent,  but  more  convenient  form: 
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where: 


Q  =  Qz  =  AcTl^^^B./U;^Cn 


P  =  P-  ^ 


The  additional  advantage  of  this  formulation  is  that  both  Q  and  H  are  in  the  stress  units 
{MPa)  instead  of  the  somewhat  cumbersome  {MPa)''.  The  directional  damage  is  defined  in 
a  manner  very  similar  to  directional  hardening,  namely 


where  u/j  are  directional  cosines  defined  in  equation  (3.14)  and  the  components  of  a  tensor 
evolve  according  to  equation 


where  q  and  M  are  material  constants.  The  initial  condition  is 

u.°(0)  =  0 

Note  that  there  are  several  problems  with  practical  application  of  directional  damage,  relia¬ 
bility  of  the  above  model  and  conducting  experiments  relevant  for  the  evolution  of  necessary 
parameters.  Even  the  extensive  experiments  presented  in  references  [11,24,25]  did  not  pro¬ 
vide  all  the  necessary  data  and,  hence,  the  damage  model  is  usually  limited  to  the  isotropic 
damage. 


3.3  Boundary  Conditions 

The  asperity  can  be  viewed  as  a  protuberance  of  a  deformable  half  space  (see  Fig.  3.12). 
It  is  subject  to  boundary  conditions  resulting  from  its  support,  contact  with  the  opposing 
surface,  adhesion  and  sliding  resistance. 

3.3.1  Support  Conditions 

If  the  asperity  is  viewed  as  the  protuberance  on  a  deformable  half  space,  the  support  condi¬ 
tions  are  defined  as  zero  displacements  at  infinity: 

lim  u  =  0 
11*11-00 
r3<0 
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In  practical  computations  we  will  usually  consider  only  a  certain  section  of  the  bulk 
material  surrounding  the  asperity.  Then  the  support  condition  will  be; 

M  =  0  on  r„ 

on  the  cut-off  boundary  Fu. 

3.3.2  Contact  Condition 

Let  the  position  of  the  rigid  flat  (see  figure  3.2)  be  defined  by; 

•  a  point  Poi^o^Voi  ^o)  which  belongs  to  the  flat,  where  Xo,yo,^o  are  its  coordinates  in 
the  initial  configuration, 

•  unit  vector  N,  normal  to  the  flat, 

•  displacement  w  of  the  flat  in  direction  N. 

Separation  of  material  point  in  the  deformable  body  from  the  flat  is  then  given  by  the 
following  formula 

s  =  x-  N  +  u-  N  —  Po-N  —  w 

where; 

X  -  initial  position  of  a  material  point, 
u  -  displacement  of  this  point. 

The  condition  that  the  asperity  can  not  penetrate  the  rigid  flat  is; 

s  >  0  on  r 

The  actual  contact  region  is  Fc  =  {x  €  F,  s(x)  =  0}.  The  difficulty  associated  with  the 
contact  condition  in  the  above  form  is  that  it  results  in  the  weak  formulation  of  the  problem 
in  the  form  of  variational  inequality,  rather  than  the  equation.  In  order  to  avoid  difficulties 
involved  in  solving  variational  inequalities,  the  contact  condition  is  usually  regularized.  In 
this  work  we  will  use  the  penalty-type  regularization  of  the  form; 

=  f^(a)  on  Fc 

where  a  =  — s  is  the  approach  (penetration)  and  is  the  value  of  traction  normal  to 
the  flat  which  defines  resistance  of  the  surface  to  penetration.  Because  we  will  be  using 
rate  formulation  in  viscoelastoplastic  analysis,  it  is  beneficial  to  introduce  a  continuously 
differentiable  penalty  function,  for  example  in  the  form  presented  graphically  in  Fig.  3.11: 
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0 

if 

a  <  0 

H  , 

if 

0  <  a  <  t 

H(a- 

if 

a  >  e 

\ 

2/ 

Here  H  is  a.  large  number  (normal  stress)  and  £  is  a  small  number  (penetration).  The 
above  penalty  function  guarantees  continuous  derivative  of  the  normal  traction  with  respect 
to  a,  which  greatly  improves  practical  performance  of  the  numerical  computations. 


3.3.3  Adhesion 

In  the  case  of  the  surface  asperity,  an  important  contribution  to  the  total  normal  force  comes 
from  adhesion  due  to  intermolecular  reactions  (see  Fig.  2.3). 

According  to  the  theory  discussed  in  our  previous  report  [94],  the  traction  boundary 
condition  due  to  adhesion  is  effective  outside  the  contact  zone: 


=  q{s)  if  s  >  0 
t%-  =  0  if  s  =  0 


where  is  the  value  of  normal  traction  OijNjNi  resulting  from  adhesion.  The  actual 
value  of  the  pressure  g  is  a  strongly  nonlinear  function  of  the  separation  s  between  surfaces 
and  according  to  the  Lennard-Jones  theory  discussed  in  reference  [65]  is  defined  by: 


q{s)  = 


8A7 

IT 


T 


(3.16) 


where  e  is  an  intermolecular  distance  (generally  e  =  0.3  —  0.5nm)  and  A7  is  the  surface 
energy  of  adhesion. 


3.3.4  Shear  Resistance 

In  addition  to  normal  stresses  on  the  contact  zone,  formation  of  metallic  bonds  provides 
shear  resistance  of  the  junction.  The  corresponding  surface  traction  can  be  expressed  in  the 
general  form: 


tp  =  tgS  +  tpT 
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where  5  is  a  vector  parallel  to  the  rigid  flat  and  T  =  N  x  S.  It  is  convenient  to  assume 
that  the  vector  S  defines  a  direction  of  prescribed  sliding  motion  of  the  surface. 

The  components  of  the  traction  vector  depend  upon  the  normal  load  (or,  equivanlently 
penetration)  and  relative  sliding  distance: 

is  =  /(a,  6s) 

h  = 

where 

a  -  penetration 

6s,  6x  -  relative  sliding  on  the  plane  in  the  S  and  T  directions  correspondingly 
6s  =  8  —  u  •  S 

h'j'  —  — 'll  •  T 

8  -  prescribed  sliding  of  the  flat  in  the  direction  of  vector  S  (see  figure  3.12) 

The  specific  form  of  functions  /  (.,.)  depends  on  the  actual  model  of  shear  resistance 
on  the  asperity  tip.  For  e.xample,  shear  resistance  governed  by  regularized  [53,59,92.93] 
Coulomb’s  friction  law  is  defined  as: 

/(a,  6)  =  /i<^^(6) 

where  is  a  coefficient  of  friction  and  (j)  is  a.  regularization  function,  defined  as: 
fi  -  friction  coefficients 

€  -  penalty  parameter  in  normal  contact  (see  section  3.3.2) 

—  1  for  b  <  Cl, 

b  6^ 

2 - h  -5-  for  Cb  <  b  <0 

6  62 

2 - r  for  0  <  b  <  Cb 

1  for  6  >  Cb 

where  Cb  is  a  small  number  (regularization  parameter).  The  form  of  the  regularization 
function  <i>  is  shown  in  figure  3.13. 

Note  that  in  practical  modeling  of  surface  asperities  the  shear  resistance  is  often  defined 
by  the  strength  of  metallic  junctions  between  asperities.  The  detailed  forms  of  corresponding 
resistance  models  will  be  studied  in  year  3  on  this  project. 

3.3.5  Initial  Conditions 

Smooth  functions  u^{x)  and  z°{x)  are  prescribed  such  that  for  a;  6  f). 
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Figure  3.13:  A  regularization  function  4> 


3.4  Variational  Formulation 

In  order  to  obtain  a  weak  formulation  of  a  boundary-value  problem  we  introduce  the  space 
of  functions 

V  =  {r  €  ,  v{x)  ^  0  as  ||®||  ^  oo} 

where  H  is  a  computational  domain,  N  is  the  dimension  of  the  physical  space  (2  or  3),  and 
is  the  Sobolev  space,  where  specific  values  of  m,p  and  q  depend  on  the  particular 
form  of  constitutive  equations. 

Multiplying  the  equilibrium  equation  (3.10)  by  a  test  function  and  integrating  over  Q.  we 
obtain  the  weak  form  of  the  rate  equilibrium  equations: 

dU  =  0  V  D  6  V' 

After  the  substitution  of  the  constitutive  equations,  application  of  the  divergence  theorem 
and  a  grouping  of  terms  the  following  variational  problem  is  obtained: 
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Find  a  displacement  rate  field  t  —*  u(x,t)  6  V  such  that 


I 


■h  /  i^N^i  "h  ■h  ixTi)vids  — 
JTc 

=  Eijkii"^lvi,jdQ  +  iiVids 


'iveV 


(3.17) 


Note  that  the  values  of  the  rates  of  normal  and  tangential  stress  on  the  contact  zone 
need  to  be  expressed  in  terms  of  displacements  using  formulas  presented  in  previous  sections 
(contact  condition,  and  sliding  resistance).  After  doing  this,  we  obtain  a  more  detailed 
formulation  of  the  problem: 


EijkiVi^jUkjdQ. 
[  E,jkiVi.ji^k}d^ 

/  ViU  ds 
Jr, 


+  ^t(5'y  +  Bfj  +  Bj^)ujds  ds  — 

+  /  Vi{G^  +  Gf  +  GJ )  ds+ 

Jr  c 

V  w  €  V'' 


(.3.18) 


where 


b;  = 

Bfi  =  U{a,bs)NiS,  +  ^{a.bs)S.S, 

Bl  =  %{aM)NiT,  +  %(aM)SiT, 

Gf  =  itiiiVj 

Gf  =  +  !{(<■.  is)  •SIS, 

Gf  =  If  (0,  h)<i’Ti 


i 


traction  from  static  boundary  conditions. 


The  rates  of  nonelastic  strains  can  be  obtained  from  the  relevant  constitutive  theory 
(see  Section  3.2.2).  For  elastic  materials  they  are  identically  equal  to  zero. 

For  elastic  materials  it  is  also  possible  and  computationally  more  efficient  to  use  a  total 
formulation  of  the  problem  which  is  as  follows. 
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Find  a  displacement  field  «(x)eV,  such  that 


/  EijkivuUkjdQ  +  /  Vi-NiNjUjds  = 

Ju  JCc  t 

Vi^Ni{w  -  Xk  ■  Nk  +  Pok  ■  Nk)ds  +  ViUds  V  r  G  K 


(3.19) 


3.5  Solution  Method  for  Elastic  Contact  Problems 

Formulation  of  the  contact  problem  is  nonlinear  even  in  the  case  of  contact  with  an  elas¬ 
tic  body  because  the  area  of  contact  depends  on  displacements  (Pc  =  r<.(u)).  Generally, 
variational  equation  (3.19)  for  u  can  be  written  in  the  following  form 

L(w)-H  =  0,  (3.20) 

where  L  stands  for  the  left-hand  side  and  R  for  the  right-hand  side. 

This  equation  is  obtained  from  (3.19)  by  requiring  that  it  be  satisfied  for  any  v  belonging  to 
V.  Because  the  condition  is  nonlinear  it  results  in  a  nonlinear  system  of  equations.  In  aim  to 
solve  the  problem  effectively,  Newton- Raphson  iteration  technique  was  used.  The  idea  of  the 
method  is  to  substitute  a  nonlinear  functional  by  its  linear  part.  Linearization  is  made  at  a 
series  of  points.  Each  point  is  the  solution  of  the  problem  linearized  at  the  previous  point.  If 
the  series  is  convergent  it  is  convergent  to  a  solution  of  the  nonlinear  problem.  For  a  simple 
case  of  a  nonlinear  equation  with  one  unknown  two  steps  of  Newton- Raphson  method  are 
shown  in  Fig.  3.14. 

Beisic  formulas  of  the  Newton- Raphson  method  are  presented  below.  Let  us  assume  that 
we  know  a  field  which  is  an  approximate  solution  of  the  equation  3.20.  Uq  =  0  can  be 
assumed.  First  two  components  of  the  Taylor  series  evolution  of  left  hand  side  of  equation 
3.20  give 


£(u")  —  R  —  gradit[jL(it")  —  R]8u  ^  0 
Assuming  that  Su  =  —  u"  we  obtain 

£(«")  -R  +  grad„L(«")(u'*+‘  -  u")  =  0 

Equation  3.20  is  a  linear  equation  for  u”"*"*.  We  use  this  equation  to  compute  a  sequence 
of  approximate  solutions  The  process  is  stopped  when  and 

\\L{u^)  —  R\\  are  small  enough. 

For  contact  with  elastic  bodies,  equation  (3.20)  has  the  following  form 

f  EijkiVijUkydQ+ f  u,-iV,wy,u”‘*'Vs  =  /  Vi-{Pok-Nk-Xk-Nk-\-w)N,ds-\- f  v,i,ds 
Jn  Jrdu’')  e  ■'  Jr^u”)  e  Jrt 
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Variations  were  evoluated  neglecting  dependence  Fc  =  rc(it).  The  above  linearized  prob¬ 
lem  is  solved  by  the  standard  FEM. 

In  order  to  provide  automatic  control  of  the  performance  of  nonlinear  procedures,  an 
expert  system-like  approach  has  been  applied.  This  application  is  based  on  our  previous 
research  on  automation  of  computational  procedures  [4],  and  employs  several  heuristic  rules 
to  monitor  and  control  the  performance  of  nonlinear  iterations.  While  in  the  original  im¬ 
plementation  discussed  in  reference  [4]  the  specialized  knowledge  engineering  software  was 
used  to  develop  the  expert  system,  in  this  project  the  essential  features  of  the  expert  system 
were  coded  in  FORTRAN  and  included  in  the  code.  The  expert  system  is  activated  at  each 
time  step  after  completing  a  prescribed  number  of  iterations  (sufficient  to  estimate  trends  in 
error  histories).  The  decisions  of  the  expert  system  are  used  to  control  the  solution  process 
and  obtain  a  converged  solution  at  minimum  cost. 


3.6  Solution  Method  for  Viscoplastic  Contact  Problems 

Formulation  of  the  problem  in  this  case  is  time-dependent. 

The  strategy  employed  in  the  solution  of  this  problem  is  as  follows:  with  the  initial 
distribution  of  stress,  temperature  and  internal  variables  specified  use  the  rate  form  of  the 
equilibrium  condition  (Eq.  (2.23))  to  obtain  the  nodal  displacement  rates.  Then  integrate 
the  constitutive  equations  forward  in  time  at  the  element  Gauss  integration  points.  With 
updated  value  of  the  stress,  temperature  and  internal  variables  at  the  new  time,  the  equilib¬ 
rium  equation  is  solved  again.  This  sequence  of  determining  the  nodal  displacement  rates, 
then  advancing  the  constitutive  equations  in  time  is  continued  until  the  desired  history  of 
the  initial  boundary-value  problem  has  been  obtained. 

Thus,  the  algorithm  proceeds  through  the  following  steps: 

1.  At  time  t,  initialize  cr,y,  Z,  for  each  element; 

2.  Calculate  e,”  =  fijicrij,Zk)  for  each  element; 

3.  Assemble  and  solve  [K](J  =  F; 

4.  Calculate  c,j  for  each  element,  e  =  [B\U\ 

5.  Calculate  d-.j  for  each  element,  a  =  [E\i  —  e"; 

6.  Calculate  Zi  for  each  element,  Zi  =  gi{crij,  Zk)', 

7.  Integrate  <7,j,  Z  forward  for  each  element  to  get  (7,j  and  Z,  at  t  -|-  At,; 

8.  If  t  -f-  At,  <  tfinai  go  to  2,  otherwise  stop. 

The  computational  method  above  has  been  presented  for  a  constant  time  step  Ats. 
Computational  experience  by  several  investigators  (see  refs.  [7,8,54])  indicates  that  a  very 
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small  time  step  can  be  required  becaus  '  nf  the  “stiff’  nature  of  the  ordinary  differential 
equations  describing  the  internal  state  variables.  To  gain  improved  efficiency  and  reliability 
a  variable  time  step  algorithm  has  been  implemented.  The  basic  idea  of  this  variable  time 
step  algorithm  is  presented  below  for  a  scalar  evolution  equation. 

The  solution  is  advanced  using  a  predictor- corrector  scheme.  The  predictor  phase  consists 
of  an  Euler  step: 

(3.21) 

The  solution  is  advanced  using  a  predictor-corrector  scheme.  The  predictor  phase  consists 
of  an  Euler  step: 

y^^+At  =  vt  +  ^tyt  (3.22) 


yt+At  ~  f  {yt+Ati^  "b 

An  error  indicator  E  [15,24]  is  then  computed  from 


(3.23) 


{yf+At  -  yt) 
^\yf+At\ 


(3.24) 


The  error  indicator  is  next  compared  with  a  preset  error  criterion  and  if  the  criterion  is 
met,  the  time  step  is  small  enough  to  proceed  to  the  corrector  stage.  Otherwise,  the  predictor 
phase  for  Eqs.  (3.22)-(3.23)  is  repeated  with  a  smaller  time  step.  For  the  viscoplastic 
evolution  equations  with  damage  modeling,  the  control  variables  used  to  calculate  the  error 
indicator  were  the  components  of  a  stress  tensor  <7,j,  internal  state  variables  Z,,  and  the 
damage  variables  w,,  with  the  maximum  of  these  selected  as  the  controlling  error. 

The  corrector  phase  is  the  modified  Newton  scheme. 


yavg  —  (jilt  +  ilt+At)  /2 


yf+At  =  yt  +  ^tyavg 

A  flowchart  depicting  the  adaptive  scheme  is  shown  in  Fig.  5.18. 
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4  Finite  Element  Analysis  of  Contact  Problems  with 
Friction 

4.1  General  Information  About  the  3D  Finite  Element  Code 

The  development  of  asperity  modeling  capabilities  was  based  on  existing,  state-of-the-art 
three-dimensional  adaptive  finite  element  kernel  code,  which  consists  of  several  separate 
modules  organized  around  the  common  data  structure  and  execution  supervisor.  Figure 
4.15  shows  a  general  structure  of  the  code  directories.  The  most  important  modules  are: 

•  an  object-based  data  structure  designed  specifically  for  the  h-p  adaptive  finite  element 
method, 

•  an  execution  supervisor  controlling  the  overall  execution  of  the  computations, 

•  pre-  and  postprocessors, 

•  an  adaptive  package, 

•  linear  equation  solvers,  and 

•  a  solver  for  a  specific  boundary  value  problem  (in  this  case,  asperity  modeling). 

Importantly,  the  elements  of  the  finite  element  data  structure,  adaptive  package,  graphical 
interfaces  and  linear  equations  solvers  were  designed  to  be  applicable  to  a  general  class  of 
problems,  and  relatively  easy  customizable  to  specific  problems  in  solid  mechanics  or  fluid 
mechanics.  Below,  selected  modules  of  the  above  kernel  are  discussed  in  more  detail. 

Object-based  Data  Structure 

A  new  state-of-the-art  data  structure  was  designed  and  implemented  in  the  kernel  to  avoid 
typical  limitations  of  traditional  finite  element  codes,  such  as: 

•  fixed  size  common  blocks  and  arrays, 

•  predefined  limits  on  problem  size, 

•  element  information  spread  throughout  memory  in  variety  of  arrays. 

The  object-based  data  structure  was  coded  in  C  computer  language,  which  allows  for 
dynamic  memory  allocation  and  more  flexible  handling  of  objects  and  structures.  Typical 
examples  of  objects  handled  by  this  data  structure  are: 

•  elements. 
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D  driver  access  routines 

Figure  4.15:  A  general  finite  element  code  structure 


47 


•  nodes, 

•  boundary  condition  data, 

•  set  of  degrees  of  freedom,  etc. 

The  major  advantages  of  object-based  handling  of  these  structures  are  listed  below: 

•  the  objects  are  created  only  when  needed, 

•  all  related  information  is  contained  in  one  structure,  and  closely  packed  in  memory, 

•  all  the  objects  are  automatically  saved/restarted, 

•  the  memory  is  reused  when  object  is  deleted,  jay  during  mesh  refinement. 

It  is  of  importance  to  note  here,  that  the  elements  in  the  above  data  structure  are 
grouped  and  colored  in  order  to  facilitate  vector  and  parallel  processing.  The  basic  idea  of 
this  vectorization  and  parallelization  is  presented  in  figure  4.16.  Importantly,  all  the  elements 
in  a  batch  are  of  the  same  type,  so  that  the  generation  of  element  stiffness  matrices  and  right- 
hand  sides  can  be  effectively  vectorized  by  putting  loop  over  elements  as  the  innermost  loop. 
On  the  other  hand,  since  the  elements  in  different  colors  have  no  common  nodes  or  sides,  the 
generation  of  element  of  element  matrices  and  assembly  for  different  colors  can  be  performed 
in  parallel. 

Adaptive  three-dimensional  finite  element  meshes 

The  finite  element  kernel  is  designed  to  handle  h-p  adaptive  finite  element  meshed  for 
three-  and  two-  dimensional  problems.  By  h-p  adaption  we  understand  a  finite  element 
technique,  wherein  the  elements  can  be  automatically  subdivided  into  smaller  elements  (h- 
refinement  /unrefinement)  and  the  polynomial  order  of  approximation  can  be  locally  in¬ 
creases  or  reduced.  A  major  advantage  of  properly  designed  h-p  mesh  is  that  it  can  achieve 
a  higher  order  of  accuracy  with  much  less  degrees  of  freedom  than  traditional  finite  element 
methods.  Moreover,  the  optimal  mesh  is  designed  automatically  by  adaptive  procedure 
driven  by  appropriate  error  estimators. 

For  three-dimensional  problems,  anisotropic  A-refinement  can  offer  a  wide  improvement  in 
computational  effort  over  more  conventional  isotropic  refinement  schemes.  This  is  primarily 
true  because  anisotropic  refinement  allows  for  selected  refinement  in  the  directions  of  interest 
only  (i.e.,  directions  of  high  error).  Thus,  anisotropic  refinement  may  greatly  reduce  the  total 
number  of  unknowns  in  many  problems,  in  turn  reducing  the  required  computational  effort. 

Isotropic  refinement  implies  that  an  element  is  identically  refined  in  each  local  direction. 
For  a  hexahedral  element,  an  isotropic  refinement  is  a  division  into  two  along  each  of  the  three 
local  directions,  which  results  in  eight  sub-elements.  In  contrast,  an  anisotropic  refinement 
of  a  hexahedral  element  is  a  division  into  two  along  a  single  local  direction,  resulting,  of 
course,  in  only  two  sub-elements.  Thus,  if  solution  phenomena  is  oriented  with  respect  to 
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Figure  4.16;  Vectorization  and  parallelizaton  for  groups  and  colors 


a  particular  local  direction,  then  anisotropic  refinement  allows  for  degrees  of  freedom  to  be 
introduced  only  in  the  direction  which  actually  reduces  the  total  error.  Isotropic  refinement, 
on  the  other  hand,  would  have  introduced  degrees  of  freedom  in  all  directions,  many  of  them 
providing  little  improvement  to  the  overall  solution.  Anisot-opic  refinement  can,  therefore, 
provide  a  higher  level  of  accuracy  than  isotropic  refinement  using  the  same  number  of  degrees 
of  freedom. 

Several  examples  of  h-adapted  meshes  in  three  dimensions  will  be  shown  in  Section  6. 
Note  that  the  mesh  refinement  introduces  several  theoretical  and  numerical  complications 
into  the  algorithm,  such  as: 

•  constrained  or  ’’hanging”  nodes  between  elements  of  different  refinement  level, 

•  propagation  of  constraints  and  possible  ’’deadlocks”  in  the  case  of  directional  refine¬ 
ments  for  complex  geometries, 

•  complications  of  unrefinement  due  to  one-to-two  approximation  rule. 

Detailed  discussion  of  these  issues  is  beyond  the  scope  of  this  report.  It  is  sufficient  to 
note,  that  before  application  of  the  above  kernel  to  asperity  modeling  all  these  difficulties 
have  been  successfully  resolved  and  the  existing  kernel  offers  operational  unique  automated 
directional  refinement  capability  for  three-dimensional  hexagonal  meshes. 

Interactive  user’s  interface  and  graphical  postprocessing 

have  been  implemented  in  the  adaptive  kernel  to  enable  user-friendly  operation  of  the  code 
and  viewing  of  three-dimensional  result.  The  interactive  graphic  interface  is  based  on  a 
window  environment,  with  a  menu-driven  selection  of  options.  A  sample  view  of  the  screen 
with  several  windows  open  is  presented  in  figure  4.17. 

The  graphic  interface  can  be  customized  for  specific  applications,  such  as  contact  and  fric¬ 
tion  modeling,  so  that  the  solution  process  and  essential  data  can  be  controlled  interactively 
by  the  user. 

The  most  important  feature  of  the  graphical  user  interface  is  a  three-dimensional  inter¬ 
active  postprocessing  capability.  For  two-dimensional  models,  such  visualization  is  rather 
trivial  as  all  of  the  computational  domain  is  always  visible  and  it  is  a  simple  matter  to  zoom 
and/or  pan  through  the  mesh  to  closely  review  the  results.  The  real  challenge  comes  from 
the  need  to  visualize  phenomena  in  three-dimensional  domains  where  most  of  the  numerical 
data  is  actually  hidden  from  the  observer  and  one  needs  to  enter  the  domain  to  view  the 
local  structure  of  the  solution. 

The  postprocessing  capability  implemented  in  the  kernel  is  capable  of  displaying  solution 
obtained  on  structured  and  unstructured  meshes,  with  both  h-refinement  and  p-enrichment 
present  in  the  mesh.  The  package  is  fully  interactive  and  operates  efficiently  on  high-end 
workstations.  The  basic  graphic  features  displayed  include: 

•  mesh  plots. 
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Figure  4.17:  Sample  screen  with  interactive  window  environment 
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•  isosurfaces  of  selected  quantities, 

•  slicing  planes  with  overlaying  isolines, 

•  deformed  configurations, 

•  three-dimensional  cursor  for  picking  pointwise  values  of  the  solution, 

All  the  above  displays  are  available  with  interactive  translation,  rotation  and  zoom  op¬ 
tions,  hidden  line  removal,  panning,  etc. 

Importantly,  the  graphics  package  is  designed  to  take  advantage  of  specialized  graphic 
hardware  and  software  available  on  many  platforms.  The  primary  platform  for  the  package 
is  the  SGI  Iris  family,  which  is  also  a  primary  platform  in  this  project.  Alternatively,  X- 
windows  graphics  is  supported,  which  is  operational  on  most  Unix  workstations. 

4.2  Formulation  of  a  Structural  Deformation  Problem  in  the  3D 
Code 

In  the  second  year  of  the  friction  project,  the  .3D  kernel  was  customized  to  solve  contact 
problems.  It  was  supplemented  with  over  5,000  instructions.  They  enable  to  run  specialized 
drivers  when  a  contact  problem  is  to  be  solved.  The  drivers  solve  the  contact  problem  using 
either  total  formulation  or  incremental  formulations  (see  section  3.4).  The  first  one  uses 
Newton-Raphson  iterative  method  and  calls  the  FEM  linear  solver  at  each  iteration  step. 
The  second  driver  uses  Euler  predictor  corrector  integration  method  with  automatic  time 
step  control  and  calls  the  FEM  linear  solver  twice  at  each  time  step. 

To  solve  solid  mechanics  problems  with  contact,  additional  customization  of  the  kernel 
FEM  code  had  to  be  implemented.  They  define,  in  a  special  format,  coefficients  of  the  volume 
and  boundary  integrals  introduced  in  previous  section.  As  regards  the  contact  condition, 
it  was  aissumed  that  contact  can  take  place  at  any  point  of  the  boundary  on  which  static 
boundary  conditions  are  applied.  Therefore,  while  the  integrals  over  this  part  of  boundary 
are  evaluated,  the  program  examines  whether  an  integral  point  is  in  contact  with  the  flat. 
If  penetration  is  greater  than  zero,  then  integrals  corresponding  to  contact  and  friction  are 
added  to  the  coefficients  of  the  stiffness  matrix  and  the  right-hand  side. 

If  a  problem  with  contact  is  to  be  solved  then  an  additional  data  file  is  read.  This  addi¬ 
tional  data  defines  initial  position  and  orientation  of  the  rigid  flat,  prescribed  displacements 
of  the  flat,  material  constants,  regularization  parameters  and  error  tolerances.  The  results 
of  each  iteration  or  time  step  are  printed  on  screen  and  to  a  disk  file. 

To  enable  automatic  generation  of  meshes,  five  additional  programs  were  prepared.  They 
generate  customized  grid  files  for; 

•  3D  axisymmetric  asperity  (cosine  hill), 

•  3D  axisymmetric  asperity  (spherical), 
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•  2D  asperity  (cosine  hill), 

•  2D  asperity  (cylindrical), 

•  2D  trapezoidal  asperity. 

The  first  two  programs  make  use  of  an  in-house  GAMMA3D  mesh  generator  .  All  of 
them  provide  generation  of  meshes  with  first  and  second  order  of  geometry  approximation. 

5  Basic  Verification  of  Numerical  Models 

To  confirm  reliability  of  the  code  and  material  models  several  numerical  tests  were  done. 
Selected  tests  are  described  in  this  section. 

The  first  test  was  carried  out  to  verify  the  total  formulation.  It  was  uniaxial  tension  of 
an  elastic  body.  It  was  performed  for  both  3D  and  2D  problems  with  both  second  and  first 
order  geometry  approximation.  The  results  of  tests  were  compared  with  analytical  solutions 
of  the  tension  problem.  The  errors  of  computed  displacements  and  stresses  were  less  then 
0.1  percent. 

The  objective  of  the  next  group  of  tests  was  to  verify  incremental  formulation  of  the 
viscoplastic  problem.  They  were  carried  out  for  alloy  B1900-|-Hf  at  temperature  871°C. 
Material  constants  as  well  as  experimental  results  for  this  material  are  given  in  reference 
[10].  For  Bodner-Partom  model  they  are  as  follows: 


Do  = 

10-*  s-i 

mi  = 

0.270  MPa-^ 

n  = 

1.03 

m2  = 

1.52  MPa-^ 

2o  = 

2400  MPa 

ri  = 

r2  =  2 

22  = 

2400  MPa 

Ai  = 

A2  =  0.0055  s 

Zi  = 

3000  MPa 

E  = 

142  GPa 

23  = 

1150  MPa 

1/  = 

0.0805 

Some  of  these  tests  are  listed  below; 

(a)  Solution  of  the  uniaxial  tension  for  an  elastic  body  by  the  incremental  formulation. 
The  results  were  the  same  as  obtained  by  the  total  formulation. 

(b)  Solution  of  the  uniaxial  tension  for  an  elasto-visco-plastic  body.  The  results  were 
compared  with  an  experiment  presented  in  reference  [55].  Certain  discrepancy  of 
results  was  observed  (see  Fig.  5.19).  This  discrepancy  was  caused  by  an  erroneous 
value  of  the  Young  modulus  given  in  reference  [55].  After  a  correction  of  Young 
modulus  {E  =  132  GPa)  the  numerical  and  experimental  results  agreed  satisfac¬ 
torily  (see  Fig.  5.20).  The  test  verified  mathematical  and  numerical  models  for 
this  simple  loading  (uniaxial  tension). 
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(c)  Next  three  tests  were  carried  out  in  order  to  examine  how  the  computer  program 
works  for  more  complicated  loading  histories.  The  results  of  those  tests  were  not 
compared  neither  with  experimental  nor  with  analytical  solutions  but  they  look 
reasonable  and  they  confirmed  reliability  of  the  computer  code.  The  tests  were 
the  following: 

•  Uniaxial  cyclic  tension  and  compression  Fig.  5.21. 

•  Uniaxial  loading  for  5,000  s  and  relaxation  for  next  5,000  s  Fig.  (5.22). 

•  Loading  for  1,000  s  and  creep  for  next  1,000  s  Fig.  5.23. 

The  above  basic  tests  verified  the  formulation  of  the  contact  problem  as  well  as  its 
applications  in  the  computer  code.  It  was  possible  to  perform  further  tests  and  comparisons 
of  results. 
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Figure  5.19:  Experimental  and  numerical  results  of  uniaxial  tension  (incorrect  value  of  Young 
modulus). 
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Figure  5.20:  Experimental  and  numerical  results  of  uniaxial  tension  of  a  specimen  with 
corrected  Young  modulus 
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Figure  5.21:  Cyclic  loading  test  of  the  viscoplastic  model  -  10  cycles. 
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6  Verification  of  Numerical  Models  of  Asperity 


Numerical  modeling  of  response  of  surface  asperities  to  contact  and  friction  loads  is  one  of 
basic  components  of  the  new  asperity-based  interface  models  developed  in  this  project.  In 
order  to  verify  correctness  of  our  finite  element  asperity  simulations,  we  performed  several 
tests  and  comparisons,  in  particular: 

•  modeling  of  elastic  sphere  in  contact  with  a  rigid  flat,  (Hertz  problem).  Analytical 
solution  is  available  for  this  problem  [47,90], 

•  numerical  modeling  and  experimental  measurements  for  custom-made  asperities  with 
strongly  pronounced  nonelastic  properties. 

Details  of  these  tests  are  presented  further  in  this  section. 


6.1  Elastic  Sphere  in  Contact  with  a  Rigid  Flat 

In  order  to  verify  the  contact  algorithm  and  the  nonlinear  solution  procedure,  a  finite  element 
solution  was  obtained  for  the  contact  of  elastic  sphere  with  a  rigid  flat.  The  finite  element 
solution  of  this  problem  was  compared  with  theoretical  solution  due  to  Hertz  [47,90].  For 
a  given  sphere  of  radius  R  and  prescribed  normal  displacements  of  the  flat  equal  iv,  the 
theoretical  predictions  of  the  contact  radius  r,  contact  area  A  and  total  load  P  are  given  by: 


r 

A 

P 


y/Rw 

wRw 


3(l-./2) 


R^wi 


The  above  problem  was  solved  numerically  using  the  following  data: 

R  =  1.0 

E  =  1000. 

1/  =  0.3 

w  =  .001,  .002,  .005,  .01,  .02 
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Figure  6.24:  Refined  mesh  for  the  Hertz  problem,  w  =  .02,  deformed  configuration.  Only 
boundary  elements  shown. 
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Figure  6.25:  Isosurfaces  of  the  vertical  displacement  for  the  Hertz  problem.  w=.02 
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Due  to  a  localized  nature  of  the  contact,  the  finite  element  mesh  was  defined  only  for  a 
section  of  the  sphere  around  the  contact  zone.  The  problem  was  solved  using  the  Newton 
procedure  combined  with  adaptive  mesh  refinement.  An  example  of  the  final  refined  mesh 
obtained  for  w=.02  is  shown  in  figure  6.24  (deformed  configuration  is  displayed  and  only 
boundary  elements  are  shown  for  clarity).  The  same  mesh  with  isosurfaces  of  vertical  dis¬ 
placement  is  presented  in  figure  6.25,  the  slicing  plane,  with  stress  cr^y,  is  shown  in  figure  6.26, 
and  the  error  indicators  projected  on  two  slicing  planes  are  displayed  in  figure  6.27. 

The  results  obtained  numerically  compare  favorably  with  numerical  predictions.  A  de¬ 
tailed  comparison  of  theoretical  and  numerical  results  are  shown  in  table  6.28,  and  the 
graphical  comparisons  of  predicted  contact  area  and  total  load  are  shown  in  figure  6.29. 
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Figure  6.27:  Distribution  of  error  indicator  for  the  Hertz  problem,  w  = 
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ELASTIC  SPHERE  IN  CONTACT  WITH  A  RIGID  FLAT 


E  = 

l.OOOOE+03 

nu=  0.3000 

R=  1.0000 

THEORY 

Displ 

Cont.  rad. 

Area 

Load 

Max.  pres. 

0.00200 

0.00500 

0.01000 

0.02000 

0.04472 

0.07071 

0.10000 

0.14142 

0.006283 

0.015708 

0.031416 

0.062832 

1.3105E-01 

5.1803E-01 

1.4652E+00 

4.1442E+00 

3.1286E+01 

4.9468E+01 

6.9958E+01 

9.8936E+01 

NUMERICAL 

Displ 

Cont.  rad. 

Area 

Load 

Max.  pres. 

0.00200 

0.00500 

0.01000 

0.02000 

0.04237 

0.07001 

0.10092 

0.14809 

0.00564 

0.0154 

0.0320 

0.0689 

1.5400E-01 

5.5400E-01 

1.5400E+00 

4.5000E+00 

4.3600E+01 

5.4000E+01 

7.9200E+01 

1.0840E+02 

Figure  6.28:  Comparison  of  theoretical  and  numerical  results  for  the  Hertz  problem 
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6.2  Experimental  Studies  of  Models  of  Asperity 

In  order  to  verify  numerical  simulation  of  nonelcistic  behavior  of  asperities,  several  experi- 
mentai  measurements  were  performed  and  then  compared  with  numerical  predictions. 

These  tests  included  simple  tension  and  compression  problems  designed  to  verify  nonelas¬ 
tic  material  constants  for  aluminum,  as  well  as  contact  tests  for  two  types  of  custom-made 
asperities.  In  this  phase  of  experiments,  which  we  refer  to  as  Phase  I,  custom  asperities  were 
chosen  in  order  to  eliminate  random  surface  factor  from  the  comparisons.  In  Phase  II,  which 
will  be  performed  in  the  next  year  of  the  project,  real  random  surfaces  will  be  considered. 

The  objective  of  this  Phase  I  experimental  study  is  to  study  the  deformation  of  contact 
surface  asperities  and  to  verify  the  analytical  prediction  by  experimentation.  In  design  of  the 
experimental  study,  it  is  initially  conceived  that  the  test  is  to  be  carried  out  under  a  small 
normal  load  (500  Lbf)  condition.  A  controlled  surface  asperities  are  machined  onto  both 
surfaces  of  an  aluminum  block.  In  the  test  arrangement,  the  aluminum  block  is  sandwiched 
between  two  hardened  steel  blocks  with  smooth  surfaces,  and  the  steel-aluminum-steel  block 
assembly  is  compressibly  loaded  to  approximately  500  Lbf.  The  deformation  of  asperities 
on  the  aluminum  block  is  monitored  during  loading.  After  completing  the  normal  loading, 
a  horizontal  load  is  slowly  applied  to  the  aluminum  block  until  the  block  begins  to  slip. 
An  estimation  of  the  coefficient  of  friction  can  thus  be  made  by  using  the  measured  normal 
and  horizontal  loads,  and  be  compared  with  that  from  the  analytical  prediction.  A  test 
apparatus  is  built  for  this  purpose,  and  tests  are  carried  out  with  this  apparatus.  After 
reviewing  the  test  results,  it  is  decided  that  the  deformation  of  surface  asperities  on  the 
aluminum  specimen  is  too  large  for  the  purpose  of  model  verification.  The  shape  of  asperity 
is  changed,  and  the  tests  are  carried  out  under  a  normal  load  of  approximately  2.000  Lbf. 
The  results  from  both  arrangements  are  reported  in  this  section. 

6.2.1  The  Test  Apparatus 

A  sketch  of  the  apparatus  is  shown  in  Fig.  6.30.  An  aluminum  (6061,  T4)  block  of  1’  x  1’  x 
0.25”  is  sandwiched  between  two  steel  blocks  of  the  same  dimension  as  shown  in  the  figure. 
The  V-shaped  grooves  of  angle  45  deg,  pitch  spacing  of  0.1  inch,  and  depth  of  0.125  inch 
are  machined  with  a  specially  designed  cutter  on  both  surfaces  of  the  aluminum  block  to 
simulate  the  surface  cisperities.  A  sketch  of  the  grooves  is  shown  in  Fig.  6.31.  The  surfaces 
of  the  steel  blocks  are  machined  smooth  and  (water)  quench  hardened  to  RC-30.  It  should 
be  mentioned  that,  due  to  the  angle  of  the  cutting  tool  and  the  pitch  spacing,  the  tips  of  the 
grooves  are  not  in  a  plane,  the  heights  of  tips  vary  alternatively  as  shown  in  the  photograph 
(Fig.  6.31)  to  be  discussed  in  a  later  section. 

A  normal  load  is  applied  to  the  specimen  assembly  through  a  mechanical  screw  jack  from 
the  bottom  of  the  apparatus.  The  normal  load  is  monitored  with  a  small  “load  transducer” 
of  capacity  500  Lbf.  The  normal  deformation  of  the  simulated  asperities  on  aluminum  block 
is  measured  with  a  “proximate  sensor”  with  an  operation  range  from  0  to  0.1  inch.  The  sand¬ 
wiched  specimen  assembly  is  first  installed  in  position,  and  an  initial  load  of  approximately 
50  Lbfs  is  applied  to  the  specimen  prior  to  “zero  adjustment”  of  the  recording  instrument 
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Figure  6.29:  Comparison  of  theoretical  and  numerical  results  for  the  Hertz  problem;  (a) 
displacement  versus  contact  area,  and  (b)  displacement  versus  contact  load 
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(an  X-Y  plotter).  A  total  normal  load  of  500  Lbs  is  then  applied  to  the  specimen  at  a  rate  of 
approximately  10  Lbf  per  minute.  After  the  normal  load  hcis  reached  500  Lbf,  a  horizontal 
load  (by  lead  blocks  and  beads)  is  slowly  applied  to  the  aluminum  block  until  the  block  starts 
to  slide.  The  horizontal  load  is  monitored  with  a  ring-shaped  load  transducer  (laboratory 
built)  as  shown  in  Fig.  6.30. 

Note  that  even  with  very  precise  calibration  of  the  test  apparatus,  certain  compliance 
or  “setting  in”  occurs  during  loading  and  pollutes  the  measurements.  This  is  especially 
true  in  case  of  high  loads  and  very  small  displacements  considered  in  this  experiment.  It  is 
a  standard  practice  in  experimental  tests  to  discard  the  initial  part  of  the  load  curve  and 
appropriately  translate  the  remaining  part.  In  this  section,  we  present  results  in  the  “raw” 
form,  but  numerical  comparisons  refer  to  corrected  graphs. 

6.2.2  Tests  Results  From  the  Above  Apparatus 

Two  tests  are  carried  out.  A  representative  force-deformation  plot  is  shown  in  Fig.  6.32.  It 
is  seen  that  there  is  an  appreciable  amount  of  plastic  deformation  of  the  surface  asperities 
on  the  aluminum  block  (neglect  the  bulk  deformation  of  both  aluminum  and  steel  blocks) 
when  the  block  assembly  is  loaded  to  500  Lbf.  A  horizontal  force  is  then  slowly  applied  to 
the  aluminum  block.  It  is  visually  observed  that  the  aluminum  block  starts  to  slip  when  the 
horizontal  load  reaches  108  Lbf.  Since  there  are  two  contact  surfaces  between  the  aluminum 
block  (with  a.sperities  on  both  surfaces  as  shown  in  Fig.  6.31)  and  steel  blocks  (with  smooth 
surfaces)  in  the  test  arrangements,  the  nominal  coefficient  of  friction  is  (108  /2)/500  =  0.108. 

6.2.3  Deformation  of  Asperities  Under  a  Larger  Normal  Load 

The  asperities  shown  in  Fig.  6.31  have  pointed  tips,  the  deformation  of  tip  is  difficult  to 
calculate.  It  is  then  decided  to  change  the  geometry  of  asperity  as  those  shown  in  Figs. 
6.31(b)  and  6.31  (c).  In  Fig.  6.31  (b),  the  asperity  is  modeled  as  a  45-deg  grooved  with 
a  truncated  tip.  In  Fig.  6.31  (c),  the  asperity  is  modeled  by  spaced  circular  rods.  The 
force-deformation  relationship  for  the  specimen  with  asperities  as  shown  in  Fig.  6.31  (b) 
and  6.31  (c)  are  measured.  The  test  arrangements  are  the  same  as  that  described  in  the 
previous  section.  Since  the  maximum  load  in  these  test  are  much  higher  than  500  Lbf.  the 
tests  are  carried  out  and  the  representative  force-deformation  curves  are  shown  in  Figs.  6.33 
and  6.34,  respectively.  A  photograph  of  the  cross-section  of  the  deformed  grooves  is  shown  if 
Fig.  6.35.  It  is  seen  that  only  the  alternate  grooves  were  deformed.  There  are  seven  grooves 
on  each  surface  of  the  aluminum  block,  therefore  only  eight  (four  on  each  surface)  of  them 
are  deformed.  Since  the  maximum  horizontal  load  required  to  move  the  aluminum  block  in 
this  exceeds  the  capacity  of  the  ring  load  cell  used  in  previous  section,  the  horizontal  pulling 
test  has  yet  to  be  carried  out  (we  need  to  build  a  new  load  cell  and  loading  frame). 

Finally,  a  test  for  the  compressive  property  of  the  aluminum  used  in  this  study  is  carried 
out.  The  stress-strain  curve  from  a  0.25”  dia  x  1”  long  aluminum  (6061.  T4)  specimen  is 
shown  in  Fig.  6.36.  Again,  the  test  is  carried  out  at  a  slow  loading  rate  (approximately  20 
Lbf  per  minute). 
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Figure  6.31:  Different  models  of  asperities  studied  experimentally. 
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Figure  6.35:  A  photograph  of  deformed  asperity 

6.3  Numerical  Simulation  of  Experimental  Measurements 

6.3.1  Viscoplastic  Uniaxial  Stress  State  Numerical  Test 

The  aim  of  this  test  was  verification  of  viscoplastic  material  constants  for  aluminum  6061 
T4.  Bodner-Partom  model  of  viscoplastic  materials  uses  14  material  constants.  However,  for 
aluminum  alloys  at  room  temperature  the  model  can  be  simplified  and  then  only  7  material 
constants  are  of  primary  importance.  Their  values  obtained  from  reference  [  ]  for  slightly 

different  heat  treatment  (T6)  are  listed  below: 

E  =  73.9  GPa 

1/  =  0.33 

Zo  =  450  MPa 

zi  =  550  MPa 

mi  =  0.12  MPa“' 

D  =  10« 

n  =  5.0 

In  order  to  verify  these  values,  the  results  obtained  with  these  constants  were  compared 
with  simple  compression  and  tension  experimental  tests.  The  test  of  uniaxial  compression 
failed  to  provide  any  useful  data  because  bifurcation  of  the  specimen  was  observed  and 
the  state  of  stresses  was  not  uniaxial.  The  adequate  verification  of  material  constants  was 
possible  using  the  results  of  tensile  test.  They  indicated  that  the  Young  modulus  was  in  fact 
smaller  and  the  yield  limit  greater  than  assumed  initially  (the  differences  reached  about  10 
percent).  In  order  to  have  a  more  accurate  model  the  Young  modulus  and  the  kinematic 
parameter  were  changed.  The  new  values  are 
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E  =  65.0  GPa 
n  =  5.8 

So  modified  material  constants  were  used  in  further  computations.  Fig.  6.38  shows  the 
plot  of  stress-strain  relation  for  initial  and  modified  material  constants  in  comparison  with 
the  experimental  results. 

6.3.2  Viscoplastic  Cylindrical  Asperity 

Aluminum  cylinder  was  used  as  a  model  of  a  "macro  asperity”.  Numerical  and  experimental 
tests  were  carried  out.  It  was  assumed  that  the  cylinder  can  be  analyzed  as  a  2D  case. 
Original  mesh  used  for  discretization  of  a  fragment  of  a  circle  is  shown  in  Fig.  6.39.  After 
the  first  computation  was  completed  the  mesh  was  automatically  refined.  It  is  shown  in  Fig. 
6.40.  The  viscoplastic  results  obtained  at  both  meshes  as  well  as  an  elastic  solution  for  the 
first  mesh  are  compared  in  Fig.  6.41.  The  conclusions  are  that: 

•  behavior  of  the  specimen  under  applied  loading  is  almost  elastic, 

•  viscoplastic  solution  is  reasonable  because  it  gives  smaller  values  of  the  contact  force, 

•  the  second  mesh  gives  results  which  can  be  treated  as  a  final  numerical  solution  (the 
diflference  between  solutions  obtained  at  coarse  and  fine  meshes  is  small.) 

The  numerical  and  experimental  solutions  are  compared  in  Fig.  6.42.  Note  that  the 
experimental  results  have  been  rescaled.  Instead  of  the  total  force  -  total  displacement 
relation  measured  in  the  experiment.  Fig.  6.42  shows  force  per  unit  length  for  upper  half¬ 
cylinder  (2  times  less  than  measured).  Moreover,  metric  units  wei’e  used. 

It  can  be  observed  that  the  numerical  results  compare  favorably  with  experimental  mea¬ 
surements.  Recall,  however,  that  the  experimental  data  were  translated,  to  correct  for 
settling  in  the  apparatus  (see  Section  6.1). 

6.3.3  Viscoplastic  Custom  Surface  Model 

Comparisons  of  results  for  a  model  truncated  V-shaped  of  asperity  are  described  in  this 
section.  Two  meshes  which  were  used  for  discretization  of  this  specimen  are  shown  in  Fig. 
6.43  and  6.44.  Numerical  results  for  both  meshes  are  shown  in  Fig.  6.46  and  6.45.  While  for 
the  cylindrical  model  the  behavior  of  the  material  was  almost  elastic  (maximum  nonlinear 
strain  was  less  than  0.4%),  in  this  case  the  influence  of  yielding  was  significant.  Maximum 
nonlinear  strain  Wcis  about  80%.  It  means  that  at  least  locally  solution  is  beyond  the  theory 
of  small  strains  which  we  use.  The  comparison  of  numerical  results  and  experimental  ones 
is  shown  in  Fig.  6.47. 

A  good  agreement  of  results  can  be  observed  for  both  models  of  asperity. 

Some  of  unevitable  sources  of  discrepancy  are  the  following 
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•  error  of  experimented  measurements, 

•  error  of  discretization  of  the  domain, 

•  error  of  modeling  of  inelastic  behavior  of  the  material, 

•  unknown  deformation  history  of  the  specimens,  prior  to  the  experiment, 

•  errors  of  time-integration  and  other  numerical  errors. 

Comparison  of  numerical  and  experimental  results  as  well  as  the  basic  numerical  tests 
allow  to  state  that  the  mathematical  model  and  the  computer  code  work  properly. 
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H  Figure  6.38:  Comparison  of  stress-strain  relation  behavior  for  uniaxial  tension  test. 
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Figure  6.39:  Original  discretization  of  the  cylinder 


Figure  6.40:  Refined  discretization  of  the  cylinder 
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Figure  6.42:  Comparison  of  numerical  and  experimental  results  for  the  cylinder 
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Figure  G.4J:  Criginal  discretization  of  the  second  specimen  (truncated  V-shaped  asperity) 


Figure  6.44:  Refined  discretization  of  the  second  specimen 
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7  A  Simple  Asperity-Based  Contact  Model 

The  results  of  elastic  solution  of  spherical  contact  problem  were  used  for  initial  tests  of 
a  complete  statistical  homogenization  procedure.  In  particular,  we  compared  our  finite 
element  based  predictions  with  a  classical  asperity-based  contact  model  due  to  Greenwood 
and  Williamson  This  is  one  of  historically  first  asperity- based  models,  derived  under 

the  following  simplifying  assumptions; 

•  the  tips  of  asperities  are  spherical, 

•  all  asperities  have  the  same  radius  of  curvature  R, 

•  asperity  height  is  a  random  variable  with  Gaussian  distribution, 

•  contact  is  elastic  and  described  by  the  theoretical  solution  of  the  Hertz  problem. 

Note  that  although  this  model  is  much  simpler  than  the  ultimate  objective  of  this  project 
(random  asperit}'  height  and  curvature,  nonelastic  deformation,  etc.),  it  provides  a  very  good 
initial  test  for  the  homogenization  procedure.  Importantly,  our  modeling  of  this  problem 
differs  from  Greenwood’s  approach  in  that: 

•  solution  of  elastic  contact  problem  has  been  obtained  from  finite  element  modeling, 
rather  than  from  Hertz  theory, 

•  expected  values  of  contact  load  and  area  have  been  calculated  using  numerical  quadra¬ 
ture,  instead  of  analytical  integration. 

Comparison  of  our  results  with  Greenwood-Williamson  theory  will  provide  a  good  mea¬ 
sure  of  errors  introduced  by  these  numerical  techniques.  Recall,  that  we  are  using  numerical 
methods  is  this  work  because  for  real  engineering  surfaces  there  exist  no  closed  analytical 
solutions  of  the  contact  problem  and  of  the  homogenization  procedure. 

The  simulation  of  the  Greenwood-Williamson  model  was  performed  for  the  following  set 
of  parameters; 

R  =  0.1414  mm 

(Tp  =  7.07106  X  lO"**  mm 

1/  =  0.3 

E  =  321.7335  kG\mm^ 
n  =  300  peaks\mm^ 

where  R  is  the  radius  of  asperity  tips,  <7p  is  the  standard  deviation  of  asperity  height  and  n  is 
the  surface  density  of  asperity  peaks.  The  results  are  presented  in  figure  7.48,  which  shows 
respective  comparisons  of  load- separation  and  load  vs. contact  area  curves.  It  can  be  seen 
that  the  difference  introduced  by  numerical  modeling  of  asperity  and  numerical  integration  of 
expectation  values  is  within  acceptable  bounds  (note  that  these  differences  could  be  further 
reduced  by  application  of  finer  meshes  for  the  asperity  contact  problem.) 
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Figure  7.48:  Comparison  of  Greenwood-Williainson  theory  with  numerical  predictions;  (a) 
load-separation  curve  and  (b)  load-area  curve. 
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8  Conclusions  and  Research  Forecast 


In  the  second  year  of  this  project  the  effort  has  been  focused  on  the  development  of  a 
finite  element  code  for  modeling  nonelastic  surface  asperities,  as  well  as  on  the  design  and 
performance  of  the  verification  experiments.  This  was  followed  by  theoretical  and  numerical 
verifications  of  the  numerical  models  of  surface  asperities.  To  model  the  deformation  of 
an  asperity,  a  generic  h/p  adaptive  finite  element  kernel  (developed  at  COMCO)  has  been 
customized  for  the  analysis  of  viscoplastic  solids  in  contact  with  a  rigid  flat,  including  elastic 
or  viscoplastic  material  properties,  contact  and  sliding  resistance. 

The  comparisons  of  numerical  predictions  of  asperity  response  with  available  analytical 
solutions  and  experimental  measurements  confirms  the  correctness  of  the  formulation  and  of 
the  finite  element  solution  procedure.  In  fact,  it  appears  that  the  actual  experiments  require 
more  caution  than  the  numerical  solutions,  because  of  calibration  difficulties  and  pollution  of 
the  results  by  the  compliance  of  the  apparatus.  These  experiences  will  be  taken  into  account 
in  the  third  year  of  this  project. 

The  effort  in  the  third  year  will  focus  on  the  actual  development  of  the  interface  mod¬ 
els  through  the  statistical  homogenization  procedure,  as  well  as  on  their  verification  and 
applications.  In  particular,  the  following  tasks  are  planned  for  the  third  year: 

1.  Evaluate  the  results  of  comparisons  performed  up-to-date  and  implement  the  necessary 
corrections  in  the  numerical  code  and  in  the  experimental  apparatus.  Design  and 
perform  experimental  measurements  for  a  random  surface  corresponding  to  a  typical 
surface  finish  (e.g.  sand  blasting). 

2.  Extend  the  finite  element  asperity  modeling  code  to  include  selected  models  of  adhesion 
and  sliding  resistance  valid  for  micro-scale  asperities  (up-to-date  the  effort  has  focused 
on  macro-scale  asperities). 

3.  Combine  the  finite  element  asperity  analysis  with  the  statistical  homogenization  soft¬ 
ware  to  develop  new  asperity-based  models  of  contact,  adhesion  and  friction.  This 
will  include  surfaces  documented  in  the  literature  as  well  as  the  surface  which  will  be 
studied  in  the  experiment. 

4.  Evaluate  the  results  of  the  above  comparisons  and  identify  possible  sources  of  discrep¬ 
ancies. 

5.  Formulate  a  theory  and  methodology  for  the  application  of  the  new  interface  models  in 
the  modeling  of  dynamic  friction  phenomena,  such  as  friction  in  presence  of  vibrations 
or  self-exited  oscillations  in  frictional  sliding. 
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